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ABSTRACT 

This  thesis  focuses  on  the  modeling  of  F-14A  component  failure  rates.  Current 
methodology  employs  the  Exponential  distribution  to  model  component  failures  and 
the  associated  Poisson  distribution  to  determine  expected  demand.  Three  other  failure 
rate  distributions  are  explored  as  alternatives:  a  WeibuU  Qight  hour  model,  a  Geometric 
sortie-dependent  model,  and  a  Mixed  sortie-Oight  hour  model.  The  expected  number  of 
component  failures  is  calculated  for  each  model  and  a  comparison  is  made  between  the 
current  model  and  these  alternatives.  The  specific  results  pertain  to  aircraft  of  this  type 
but  the  concepts  employed  can  be  applied  to  other  aircraft  as  well. 

The  Geometric  model  provided  a  better  fit  for  components  which  were  not 
operated  continuously,  and  the  WeibuU  performed  better  when  the  components  were 
operated  continuously.  Overall,  the  Exponential  was  the  least  elTective  model  for  the 
nine  components  studied. 


THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may 
not  have  been  exercised  for  all  cases  of  interest.  While  ever\'  efTort  has  been  made, 
within  the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and 
logic  errors,  they  cannot  be  considered  validated.  Any  application  of  these  programs 
without  additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION 

A.       BACKGROUND 

1.  A  Carrier  Naxy 

The  carrier  is  the  nucleus  of  today's  U.S.  Navy,  but  its  ability  to  project  power 
is  a  byproduct  of  the  synergistic  relationship  it  shares  with  its  air  wing.  As  Peter 
Garrison  stated  in  CV:  Carrier  Aviation,  "by  adopting  submarine-launched  ballistic 
missiles  and  ship-based  tactical  aircraft  the  navy  has  co-opted  some  of  the  tools  of  the 
army  and  air  force  and  adapted  them  to  a  marine  environment.  Ships  in  themselves  are 
no  longer  powerful,  they  take  their  power  from  the  airplanes  or  long  range  missiles  that 
they  carr\'.  An  adversary'  setting  out  to  sink  a  ship  does  not  send  another  ship  against 
it,  but  instead  sends  an  airplane  or  a  missile"  [Ref  1].  (italics  mine).  The  fleet  has 
become  dependent  upon  aircraft  to  take  the  battle  over  the  horizon,  and  to  extend  the 
battle  group's  area  of  influence.  Without  downplaying  the  importance  of  the  carrier,  it 
is  just  another  large  ship  without  its  contingent  of  aircraft.  Thus,  the  carrier's  real 
importance  to  the  fleet  rests  on  the  abiUty  of  its  personnel  to  support  and  maintain  the 
air  wing  at  an  optimal  level  of  readiness. 

2.  Aircraft  Readiness 

Ideally,  the  fleet  would  like  each  aircraft  to  be  fully  mission  capable  (FMC); 
that  is,  the  material  condition  of  the  aircraft  should  satisfy  minimum  requirements  to 
perform  all  of  its  missions.  An  aircraft's  mission  status  (e.g..  partially  mission  capable 
(PMC),  not  mission  capable  (NMC).  etc.)  is  determined  using  the  Mission  Essential 
Subsystems  Matrices  published  in  OPNAVINST  5442.4H. 

The  ability  to  ensure  a  high  percentage  of  FMC  aircraft  is  constrained  by 
several  factors.  First,  the  Xa\7  has  not  been  able  to.  provide  shipboard  intermediate 
maintenance  activities  (IMA)  with  sufficient  numbers  of  qualified  personnel,  especially 
in  critical  technical  billets.  In  part,  this  results  from  the  allocation  of  manpower  based 
on  squadron  size  (i.e..  the  smaller  squadrons  have  fewer  personnel  to  maintain  their 
aircraft).  The  problem  associated  with  the  manpower  shortage  is  compounded  by  a 
need  for  specialization  in  troubleshooting  and  repairing  the  different  aircraft  types, 
models,  and  series  deployed  on  the  carrier.  Second,  the  IMA's  capability  to  repair 
components  as  they  fail  is  constrained  by  limitations  of  test  equipment,  maintenance 
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facility  size,  availability  of  repair  parts,  and  the  total  demand  within  the  repair  pipeline. 
Finally,  those  components  which  cannot  be  fixed  at  the  intermediate  level  must  be 
resupplied  by  outside  sources. 

The  air  wing  must  be  able  to  accurately  identify  and  communicate  its 
requirements  for  component  replacements  and  resupply  facilities  should  be  capable  of 
expeditiously  responding  to  the  those  needs.  During  periods  of  high  tempo  flight 
operations,  as  would  be  experienced  during  wartime,  the  factors  listed  above  become 
more  critical.  The  ability  to  maintain  high  levels  of  readiness  is  hampered  by  an 
increase  in  the  number  of  components  in  the  repair  pipeline,  manpower  and  facility 
limitations,  and  a  reduced  or  nonexistent  abiUty  to  resupply  vital  components. 
Potential  problems  may  be  partially  rectified  by  improving  the  methodologies  currently 
used  to  attain  that  "best"  possible  mix  of  replacement  parts  located  on  the  aircraft 
carrier.   [Ref.  2:  pp.  3.9] 

B.       AVIATION  CONSOLIDATED  ALLOWANCES 

1.  Definition 

The  Aviation  Consolidated  Allowance  List  (AVCAL)  is  supposed  to  be  a 
"best"  mix  invehtor\'.  An  AVCAL  is  defined  in  OPNAVINST  4790. 2C  as  a  list  of 
aeronautical  material  tailored  to  each  individual  ship  to  support  assigned  or  embarked 
aircraft  flight  operations.  This  list  of  parts  is  prepared  by  the  Navy's  Aviation  Supply 
Ofiice  (ASO)  under  the  direction  of  the  aircraft  type  commander,  and  determines  spare 
part  allowances  based  on  expected  demands. 

Ideally,  it  would  be  advantageous  to  stock,  at  least  one  spare  for  each  part 
that  could  possibly  fail  during  deployment  regardless  of  how  unlikely  that  occurrence 
might  be.  Unfortunately,  inventory  sizes  are  constrained  by  both  budget 
considerations  and  storage  limitations.  The  latter  is  the  result  of  a  diverse  collection  of 
aircraft  types,  manufactured  by  different  companies  with  each  aircraft  requiring  its  own 
specific  replacement  parts,  coupled  with  the  inherent  aversion  to  building  aircraft 
carriers  even  larger  than  they  already  are  simply  to  carry  more  parts.  Table  1  lists  the 
normal  complement  of  aircraft  aboard  a  carrier,  and  is  provided  to  illustrate  the 
differences  in  aircraft  types,  missions,  and  spare  part  requirements.   [Ref  2:  p. 5]. 

2.  Need  for  Improvement 

Determination  of  the  proper  AVCAL  allowance  for  each  component  is  not  an 
easy  task,  and  is  definitely  an  area  needing  improvement.  A  recent  study  conducted  by 
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TABLE  1 

NORMAL  COMPLEMENT  OF  CARRIER  AIRCRAFT 

AIRCR.AFT 

MISSION 

SQUADRONS 

AIRCR.'XFT 

F-14A 

Fighter 

2 

24 

A-7E 

Light  Attack 

2 

24 

A-6E 

Attack 

10 

KA-6D 

Tanker 

4 

S-3A 

ASW 

10 

E-2C 

SURV 

4 

EA-6B 

EW 

4 

SH-3 

ASW 

6 

TOTAL 

81 

the  Center  for  Naval  Analyses  reported  that  "statistics  for  combat  aircraft  operating 
during  the  period  July  1982  through  June  1983  show  aircraft  operating  at  less  than  60 
percent  of  their  anticipated  wartime  rate.  These  aircraft  were  not  FVIC  due  to  supply 
between  18-24  percent  of  the  time.  These  readiness  rates  would  drop  significantly  if  the 
aircraft  flew  at  a  wartime  rate  and  AVCALs  continued  to  be  constructed  as  they  were 
in  this  period"  [Ref.  3:  p.  1-1].  This  statement  strengthens  the  argument  that  current 
methods  of  AVCAL  formulation  are,  at  best,  adequate  for  peacetime  operations,  but 
could  result  in  disastrously  low  readiness  rates  when  exposed  to  the  additional 
requirements  associated  with  wartime  utilization.  Two  approaches  to  detennining  a 
better  AVCAL  would  be  the  introduction  of  new  methods  of  modelling  requirements  or 
the  improvement  of  weak  areas  within  the  current  model. 
3.  AVCAL  Models 

Two  categories  oC  models  used  to  compute  the  AVCAL  are  readiness-based 
and  demand-based  inventor>'  systems.  Readiness-based  systems  determine  spare 
allowances  that  achieve  a  predetermined  level  of  readiness  while  minimizing  the  cost  of 
the  AVCAL.  Examples  of  readiness-based  models  are  Rand's  DYNAMETRIC  model 
and  the  Center  for  Naval  Analyses'   Ml  VIE   model  [Ref  3:  p.  2-8].     Demand-based 


14 


(supply  effectiveness)  models  determine  component  allowances  based  on  expected 
demand  and  the  capabilities  of  the  repair  facility  and  resupply  pipeline.  The  ASO 
Manual  model  is  an  example  of  a  demand-based  method,  and  is  the  model  currently 
used  by  the  Navy  to  compute  AVCAL  spare  allowances. 

C.       ASO  MANUAL  MODEL 

The  ASO  Manual  Model  uses  historic  data  to  forecast  demand  rates  and 
component  repair  times.  The  model  determines  the  type  and  number  of  components 
included  in  the  AVCAL  allowance.  For  a  detailed  explanation  of  the  rules  governing 
component  selection,  the  reader  should  refer  to  either  Aviation  Parts  Allowance  Policy 
by  Peter  Evanovich  or  A  Retail  Inventory  Model  for  lYaval  Aviation  Repairable  Items  by 
Mark  L.  Mitchell.   [Refs.  3,4]. 

1.  Attrition  Only  Items 

Attrition  only  items  are  repairable  components  that  cannot  be  repaired  at  the 
intermediate  maintenance  activity  aboard  the  carrier.  The  expected  demand  for  items 
of  this  type  is  computed  using  equation  I.l. 

NB'  =  NBx(FH'/FH)  (eqn  1.1) 

where, 

NB'  =  number  of  expected  attritions  for  a  specific  flight  program 

NB    ='  number  of  attritions  in  the  ASO  data  base  ■ 

FH'  =  specified  flight  program 

FH    =  total  flight  hours  represented  by  the  ASO  data  base 

2.  Items  that  can  be  Repaired  Locally 

The  allowance  for  items  that  can  be  repaired  locally  is  composed  of  two  parts. 
An  allowance  for  repairable  components  that  are  beyond  the  capabilities  of  the  local 
maintenance  activity  to  repair  them  (BCM),  and  an  allowance  for  components  that  are 
queued  at  various  places  in  the  local  repair  cycle. 

a.   BCM  Allowance 

Units  are  considered  BCM  if  the  maintenance  activity  cannot  repair  the 
item  because  of  a  lack  of  tools,  equipment,  parts  or  technical  expertise.  Equation  1.1  is 
used  to  forecast  the  number  of  BCM  failures  expected  to  occur  tor  a  specific  flight 
program. 
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b.   Repair  Pipeline  Allowance 

The  second  input  to  the  local  repairable  allowance  provides  units  to 
compensate  for  the  number  of  failed  components  undergoing  repair  at  the  IMA.  This 
allowance,  called  the  Local  Repair  Cycle  Allowance  (LRCA),  is  set  at  a  quantity  which 
should  ensure  that  demand  does  not  exceed  available  supply  at  least  ninety  percent  of 
the  time.  This  safety  level  is  dependent  upon  an  accurate  estimate  of  the  number  of 
components  repaired  by  the  shipboard  intermediate  maintenance  activity.  Such  a 
forecast  is  expressed  in  equation  1.2. 

NR'  =  NRx{FH7FH)  (eqn  1.2) 

where, 

NR'  =  number  of  expected  repairs  for  a  specific  flight  program 
NR    =  number  of  repairs  represented  in  the  ASO  data  base 
and, 
.  FH'  and  FH  remain  as  defined  above 

D.       FAILURE  RATE  DISTRIBUTION 

Combining  equations  1.1  and  1.2  provides  the  expected  number  of  total  failures, 
both  repairables  and  attritables,  for  the  specified  flight  program.^ 

..     '.  NR'  +  NB'  =  (NR  +  NB)x(FH';FH)  (eqn  1.4) 

Substituting  NF'  for  NR'  +  NB'  and  NF  for  NR  +  NB.  equation  1.4  can  be  expressed 
as: 

NF'  =  NFX(FH'/FH)  (eqn  1.5) 

Equation  1.6  results  from  rearranging  terms  on  the  right  hand  side  of  equation  1.5. 

NF'  =  (NF;FH)xFH'  (eqn  1.6) 


^For  "attrition  only  items"  NR'  equals  zero. 
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The  total  number  of  expected  failures  NF'  equals  the  product  of  some  specific 
future  flight  program  (FH')  and  the  quantity  NF/FH.  The  latter  should  be  easily 
recognized  as  the  maximum  likelihood  estimate  (based  on  the  historical  data  base)  for 
the  failure  rate  lambda.  Thus,  NF'  is  the  expected  value  of  a  Poisson  distribution  with 
the  failure  rate  parameter  X  and  interval  of  observation  {0,FH'].  The  use  of  the  Poisson 
distribution  implies  that  the  underlying  distribution  for  the  time  between  failures  is 
exponential.  Certain  properties  of  the  exponential  distribution  simplify  the 
computational  problem,  but  what  price  is  paid  for  this  simplicity?  The  remainder  of 
this  paper  will  suggest  alternative  distributions  for  describing  the  time  between 
component  failures  and  estimating  the  expected  number  of  failures  for  a  given  flight 
program. 

E.       SYNOPSIS 

Chapter  II  describes  the  methods  and  assumptions  used  to  extract  component 
failure  information  from  a  data  base  containing  six  squadrons  of  F-14A  aircraft  flight 
and  maintenance  records.  Chapters  111  and  IV  discuss  current  and  possible  methods 
of  modeling  aircraft  failure  rates,  and  techniques  used  to  derive  maximum  likelihood 
estimates  of  model  parameters.  Flight  hour  dependent  models  are  presented  in  Chapter 
III  and  include  the  traditional  Exponential  model  and  a  WeibuU  model.  Chapter  IV 
discusses  a  sortie-dependent  Geometric  model,  and  the  sortie-flight  hour  Mixed  model 
proposed  in  Distinguishing  the  Effects  on  Failures  of  Changes  in  Sortie  Rate  and  Sortie 
.  Leng'th  by  Robert  A.  Levy  [Ref  5].  The  four  different  models  will  be  compared  in 
Chapter  V  using  both  graphical  and  quantitative  measures  of  fit.  Finally,  Chapter  VI 
will  discuss  the  advantages  and  disadvantages  of  each  model,  along  with  the 
conclusions  and  recommendations. 
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II.  DATA  BASE 

A.  ORIGIN  OF  DATA 

The  initial  data  base  was  obtained  from  the  'Savy  Maintenance  Support  Office 
(NMSO),  Mechanicsburg,  Pennsylvania  via  the  Center  for  Naval  Analyses,  Alexandria, 
Virginia.  The  data  set  included  the  historical  flight  and  maintenance  records  for  two 
operational  deployments  of  F-14A  aircraft.  The  most  recent  deployment  included  two 
squadrons,  VF-1I4  and  VF-213,  which  operated  ofTthe  USS  ENTERPRISE  (CVN-65) 
during  the  period  of  March  1984  through  December  1984.  This  deployment  became  the 
"model"  set,  with  its  flight  and  maintenance  records  providing  the  parameter  estimates, 
for  the  study. 

The  second,  earlier  deployment  consisted  of  four  squadrons,  VF-11  VF-31  and 
VF-14  VF-32,  which  operated  olT  the  USS  KENNEDY  {CV-67)  and  USS 
INDEPENDENCE  (CV-62),  respectively,  during  the  period  of  November  1983  through 
March  1984.    In  this  study  these  four  squadrons  were  treated  as  the  "validation"  set. 

Using  the  parameter  estmiates  from  the  model  set,  four  specific  failure  rate 
distributions  were  used  to  calculate  the  expected  number  of  component  failures  which 
should  occur  during  a  deployment  similar  to  the  one  flown  by  the  four  squadrons 
aboard  the  USS  KENNEDY  and  USS  INDEPENDENCE.  The  actual  number  of 
component  failures  observed  during  the  deployment  of  these  squadrons  was  then 
compared  with  the  expected  number  generated  using  the  parameters  estimated  from  the 
"model"  set  deployment. 

B.  MAINTENANCE  DATA  SYSTEM  (MDS) 
1.  Purpose 

The  maintenance  data  system  (MDS)  was  developed  to  improve  the 
management  capabilities  of  the  Navv-'s  Vlaintenance  and  Material  Management  (3-M) 
system.  It  was  designed  as  a  management  information  system  which  would  provide  the 
statistical  data  needed  to  analyze: 

•  Equipment  maintainability  and  reliability 

•  Equipment     configuration    to     include    alterations     and    technical    directive 
compliance 

•  Equipment  mission  capability  and  utilization  rates 

•  Material  usaee 
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•  Maintenance  and  material  procurement  times 

•  Weapon  system  and  material  costing 


To  meet  the  data  requirements  listed  above,  NMSO  divides  the  data  collection 
efibrt  into  four  specific  areas:  material  reporting  (MR),  subsystem  capability  impact 
reporting  (SCIR),  maintenance  data  reporting  (MDR),  and  utilization  reporting.  The 
material  reporting  system  documents  those  materials  used  in  support  of  the 
maintenance  efTorts.  The  subsystem  capability  impact  report  provides  the  higher  level 
commanders  with  information  pertaining  to  their  subordinate  command's  ability  to 
conduct  their  mission.  The  last  two  areas  of  data  collection,  maintenance  data 
reporting  and  utilization  reporting,  provided  the  data  necessary  to  determine  the  times 
between  component  failures. 

2.  Maintenance  Data  Reporting  (MDR) 

The  documents  which  support  the  maintenance  data  reporting  system  include 
support  action  forms  (SAF),  metrology  equipment  recall  (VIETER)  cards,  and  the 
Visual  Information  Display  System/Maintenance  Action  Forms  (VIDS/MAF).  Support 
action  forms  document  maintenance  actions  that  did  not  require  any  corrective  action. 
such  as,  aircraft  servicing,  aircraft  handling,  gind  the  preliminary  look  phase  of  an 
inspection.  VIETER  cards  identify  the  maintenance  time  spent  calibrating,  and  repairing 
the  test  and  measuring  systems.  The  last  form.  VIDS/VIAF  (OPNAV  4790/60), 
documents  inspections,  technical  directive  compliance,  and  repair  actions.  It  identifies 
the  component  failures  by  work,  unit  code,  time  of  failure  (Julian  date  and.  when 
discovered  code),  aircraft  experiencing  the  failure  (bureau  number),  cause  of  failure 
(malfunction  description  code),  and  type  of  maintenance  required  (transaction  code  and 
action  taken  code). 

3.  Utilization  Reporting 

The  flight  records  (utilization  rates)  were  compiled  from  information  recorded 
on  the  Aircraft  Yellow  Sheet  (OPNAV  3760/2B).  The  Yellow  Sheet  is  prepared  by  air 
crew  personnel  at  the  end  of  each  mission  and  contains  information  identifying  the 
flight  according  to  aircraft  flown,  mission  type,  Julian  date,  flight  duration,  etc. 

C.       SELECTED  WORK  UNIT  CODES 

A  work  unit  code  (WUC)  is  an  alphanumeric  code  identifying  an  item  on  which 
work  is  being  performed.  Systems  are  identified  with  two  digit  codes  and  components 
with  Ave  digit  codes.  Seven  digit  codes  are  used  to  further  breakdown  components  into 
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lower  levels  of  subassembly.  The  components  selected  for  this  study  are  listed  in 
Table  2.  Each  of  the  nine  components  chosen  had  at  least  thirty  failures  during  the 
deployment  cycle. 


TABLE  2 

COMPONENTS  WITH  ASSOCIATED  WORK  UNIT  CODES 

COMPONENT 

WUC 

CN1263;ASN92(V)  INERTIAL  MEASURING  UNIT 

734H100 

ECU74/A  SIGNAL  DATA  CONVERTER 

46X1600 

MXU611/A  FUEL  TANK  RELEASE  MECHANISM 

4622100 

CP1166;A  AIR  DATA  COMPUTER 

56X2500 

AN/ARC159/159(V)-1  UHF  TR.-\NSMltTER 

632Z100 

T1224/AWG9  R.-XDAR  TR.ANSMITTER 

74A1500 

IP1185/AWG9  DETAIL  DATA  DISPLAY 

74A5M00 

CP124S/ASW43  AIRCR.-XFT  ROLL  COMPUTER 

5772200 

IP1027/AVA12  ANALOG  DISPLAY  INDICATOR 

6918100 

D.       DATA  REDUCTION 

The  intent  of  the  data  reduction  process  is  to  eliminate  those  maintenance 
actions  which  do  not  identify  a  component  failure  (e.g.,  aircraft  servicing, 
troubleshooting,  calibrations,  or  other  routine  maintenance  actions).  The  maintenance 
codes  not  indicative  of  a  component/item  failure  are  found  within  the  series  of 
malfunction  description,  transaction,  and  action  taken  codes.   [Ref  6:  pp.  17-19] 

1.  Malfunction  Description  Codes 

The  three  digit  numeric  malfunction  description  codes  describe,  as  accurately 
as  possible,  the  probable  cause  of  a  system; component  failure.  Vlaifunction 
description  codes  used  when  no  failure  was  detected  are  displayed  in  Table  3. 
Maintenance  records  which  contained  these  codes  were  eliminated  from  the  final  data 
set.  Conditional  malfunction  description  codes  identify  component  failures  caused  by 
factors  other  than  natural  wear-out.   Examples  of  this  type  of  failure  include:  failed  due 
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to  weather,  foreign  object  damage  (FOD),  improper  or  faulty  maimenance,  etc.  To 
reduce  the  complexity  of  parameter  estimation,  conditional  failures  were  treated  as  if 
the  failure  had  occurred  through  normal  usage.  The  components  studied  had  a  low 
percentage  of  conditional  failures.  In  fact,  less  than  three  percent  of  all  observed 
failures  were  conditional. 


TABLE  3 
EXCLUDED  MALFUNCTION  DESCRIPTION  CODES 

CODE      DESCRIPTION 

799  .   No  defect 

800  No  defect  -  removed  and/or  reinstalled  to  facilitate' maintenance 

801  No  defect  -  removed  for  modification   - 

803  No  defect  -  removed  for  time  chanse 

804  No  defect  -  removed  for  sclieduled "maintenance 

805  No  defect  -  removed  for  pool  stock 

806  No  defect  -  removed  as  part  of  a  matched  system 

807  No  defect  -  removal  directed  bv  hiaher  authority 
811  No  defect  -  removed  during  troubleshooting 


2.  Transaction  Codes 

Transaction  codes  are  two  digit  numeric  codes  used  to  identify  the  type  of 
data  reported  by  the  maintenance  activity.  Transaction  codes  used  in  the  maintenance 
reporting  system  are  listed  in  Appendix  A.  Table  4  lists  those  transaction  codes 
indicative  of  maintenance  actions  that  required  removal,  installation,  and/or  repair  of  a 
defective  component.  Maintenance  records  associated  with  these  codes  were  included  in 
the  final  data  set.  All  other  maintenance  records  were  deleted. 

3.  Action  Taken  Codes 

Action  taken  codes  are  one  digit  alphanumeric  codes  which  describe  the 
maintenance  action  performed  on  the  system,  component  or  item.  A  complete  list  of 
action  taken  codes  is  shown  in  Appendix  B.  Action  taken  codes  associated  with  the 
removal,  installation,  and/ or  repair  of  a  defective  component  were  included  in  the  final 
data  set  and  are  shown  in  Table  5. 
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TABLE  4 

TRANSACTION  CODES  USED  IN  STUDY 

CODE 

DESCRIPTION 

12 

On-equipment  work,  includine  engines,  involving  non- 
repairable  components/items  documented  as  failed  parts. 

21 

Will  be  used  when  a  repairable  component  is  removed, 
excluding  engine  components,  for  processing  at  an  I  or  D 
level  mamtenance  activity  (This  code  is  used  when 
only  the  removal  must  be  documented  and  a  replacement 
is  not  required). 

23 

Removal  and  replacement  of  a  defective  or  suspected 
defective  component  from  an  end  item,  excluding  engines 
at  the  O-level.  Additionallv  this  code  is  used  for^the^ 
removal  and  replacement  of  a  complete  engine  assembly 
for  a  defect,  suspected  defect,  or  scneduled^maintenance 
requirement. 

24 

Will  be  used  when  a  repairable  engine  component 
is  removed  for  orocessing  at  an  I  or  D  level  activity 
(This  code  is  used  when  onlv  the  removal  must  be' 
documented  and  a  replacement  is  not  required). 

25 

Removal  and  replacement  of  a  defective  or  suspected 
defective  repairaole  component  from  an  'engine. 

E.       FINAL  DATA  SET 

1.  Problems 

Several  problems  were  encountered  during  the  data  reduction  process  that 
require  discussion.  First,  it  was  not  always  clear  which  flight  evolution  caused  a 
component  to  fail.  The  cause  for  this  confusion  is  the  poor  record  keeping  associated 
with  the  actual  time  of  component  failures.  For  example,  the  maintenance  data 
reporting  system  does  not  record  the  actual  time  of  component  failure.  Instead,  the 
data  is  arranged  according  to  the  time  the  maintenance  activity  received  the  component 
for  repair.  In  most,  but  not  all  cases,  the  component  was  received  by  the  maintenance 
activity  immediately  following  the  flight  causing  the  component  failure. 

Second,  several  records  indicated  multiple  failures  of  a  single  component.  For 
example,  a  single  component  may  experience  several  unrelated  malfunctions  during  a 
single  flight  evolution.  To  expedite  the  repair  and,  or  replacement  of  the  failed 
component,  the  maintenance  reporting  system  requires  the  documentation  of  each 
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TABLE  5 

ACTION  TAKEN  CODES 

CODE 

DESCRIPTION 

1 

BCM  -  maintenance  activity  not  authorized  to  conduct 

repair 

2 

BCM  -  Lack  of  equipment,  tools,  or  facilities 

3 

BCM  -  Lack  of  technical  skills  necessary  to  complete 

repair 

4 

BCM  -  Lack  of  parts 

5 

BCVI  -  Fails  check  and  test,  and  maintenance  is  allowed 
to  conduct  check  and  test  only 

6 

BCVI  -  Lack  of  technical  data 

7 

BCM  -  Beyond  the  authorized  repair  depth 

8 

BCM  -  Administrative 

9 

BCM  -  Condemned,  repair  not  feasible 

C 

This  code  is  used  when  a  repairable  iiem  of  material 
identified  by  a  work  unit  code  is  repaired. 

F 

Failure  of  components/items  undergoing  test. 

P 

This  code  is  entered  when  an  item  of  material  is 
removed  and  only  the  removal  is  to  be  accounted  for. 

Q 

This  code  is  entered  when  an  item  is  installed  and 
only  the  installation  action  is  to  be  accounted  lor. 

R 

This  code  is  entered  when  an  item  of  material  is 
removed  due  to  suspected  malfunction  and  the  same  or  a 

like  item  is  reinstalled. 

*  BCM:  Action  taken  codes  labelled  1-9  are  used  for  repairable 
items  which  have  been  adniinistrativelv  and;  or  technicallv  screened 
and  lound  to  be  beyond  the  capability  ol  the  maintenance  activity. 

malfunction  (i.e..  multiple  failures  of  a  single  component).  Previous  studies  have 
treated  this  problem  as  independent  single  component  failures  with  identical 
interoccurrence  times.  The  alternative  is  to  treat  the  multiple  failure  as  a  single 
incident. 

Third,  the  underlying  assumption  of  the  exponential  distribution  is  that 
components  fail  as  a  function  of  flight  hours  flown.  This  would  imply  that  a 
component  is  operated  continuously,  which  may  or  may  not  be  the  case. 
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Fourth,  it  was  difficult  to  determine  if  a  component  was  replaced  immediately 
following  the  flight  evolution  causing  the  failure.  Additionally,  there  was  some  concern 
pertaining  to  periods  of  flight  operations  during  which  a  component  might  be  missing. 

Fifth,  the  condition  of  the  component  (i.e.,  new  or  used)  at  the  start  of 
deployment  was  not  always  known.    A  related  problem  involves  the  use  of  the  same 
failure  rate  distribution  for  both  repaired  and  new  components. 
2.  Simplifying  Assumptions 

Before  proceeding  with  any  data  analysis,  it  is  necessary  to  make  assumptions 
which  can  simplify  or  eliminate  the  problems  discussed  in  the  previous  section.  These 
assumptions  are  presented  in  the  same  order  as  the  problem  areas  they  are  intended  to 
simplify. 

First,  all  failures  are  associated  with  the  sortie  during  which  they  occurred  for 
When  Discovered  Codes  A,  B.  C.  and  D.  For  all  other  codes,  it  is  assumed  that  the 
training  and  experience  of  the  flight  and  ground  crews  ensure  the  discover\'  of  the 
failure  at  the  earliest-  possible  opportunity  (i.e.,  the  failure  occurred  during  the  flight 
prior  to  discover^').  A  complete  listing  of  When  Discovered  Codes  is  found  in 
Appendix  C. 

Second,  all  multiple  failures  of  a  single  component  are  treated  as  a  single 
component  failure.  It  is  believed  that  multiple  failures  would  affect  the  turnaround 
time  required  to  repair  the  component  and  may  cause  delays  within  the  repair  pipeline. 
However,  this  should  not  affect  the  rate  at  which  the  component  fails. 

Third,  it  is  assumed  that  components  operate  continuously,  at  least  in  a 
standby  mode.  If  this  is  not  the  case,  the  sortie-dependent  models  should  provide  a 
better  estimate  of  expected  failures. 

Fourth,  the  study  assumes  that  components  are  replaced  inunediately  and  that 
gaps  in  coverage  are  not  allowed.  In  fact,  it  is  extremely  unlikely  that  an  aircraft 
would  operate  without  any  of  the  components  listed  in  this  thesis. 

Fifth,  all  components  are  considered  new  at  the  start  of  the  deployment.  For 
the  Exponential  and  Geometric  models  this  assumption  is  not  necessan.'  due  to  the 
"memor\'less  property"  they  exhibit.  Additionally,  the  impact  on  the  WeibuU 
distribution  should  be  minimal  if  the  number  of  observed  failures  is  large. 

Finally,  to  reduce  the  complexity  of  the  problem,  it  is  necessar\'  to  assume 
that  repaired  components  exhibit  the  same  failure  characteristics  as  new  components 
(i.e.,  have  the  same  failure  rate). 
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3.  Format 

a.  Parameter  Estimation 

To  obtain  the  best  possible  estimate  of  the  time  between  failures,  the 
observations  within  the  initial  data  base  were  sorted  by  aircraft  bureau  number,  then 
Julian  date  and  the  event  time.  Failure  times  were  measured  between  consecutive 
component  failures  on  a  specific  aircraft.  These  measurements  were  never  made 
between  component  failures  on  two  different  aircraft.  The  information  included  in  the 
final  data  set  is  the  aggregate  totals  and  includes  a  censoring  indicator,  number  of 
sorties  flown,  flight  time  without  failure,  and  flight  time  with  failure. 

A  censor  indicator  equal  to  zero  indicates  that  a  component  failure  was 
observed.  The  data  entries  associated  with  this  type  of  observation  include  the  number 
of  flights  flown  inclusive  of  the  event  causing  the  failure,  total  flight  time  flown  since 
last  failure  (does  not  include  flight  time  of  event  causing  failure),  and  the  flight  time  of 
the  event  causing  the  failure. 

A  censor  indicator  equal  to  one  represents  an  event  which  is  right  hand 
censored.  Right  hand  censoring  refers  to  an  observation  of  a  component  which  does 
not  fail  prior  to  the  termination  of  the  deployment  (i.e.,  the  data  set  documents  a 
components  survival  time,  but  provides  no  information  pertaining  to  the  actual  time  of 
component  failure).  The  data  entries  associated  with  this  type  of  observation  are  the 
total  number  of  flights  flown  and  the  total  flight  time  since  the  last  component  failure. 

b.  Model  Comparison 

The  intent  of  this  thesis  is  to  determine  each  model's  ability  to  predict 
failures  for  the  components  listed,  and  compare  those  predictions  with  the  actual 
component  failures  observed  during  the  deployment  of  the  squadrons  aboard  the 
USS  KENNEDY  and  L'SS  INDEPENDENCE.  To  accurately  assess  these  capabilities, 
it  was  necessary'  to  identify  component  failures  as  they  occurred  during  the  deployment. 
Thus,  the  data  set  used  to  provide  the  graphical  and  quantitative  comparison  of  the 
predicted  vs  actual  component  failures  was  different  than  the  data  set  used  to  estimate 
the  model  parameters.  The  latter  data  set  was  sorted  by  Julian  date  and  event  time 
only.  The  accumulated  number  of  sorties  flown  and  flight  hour  totals  were  recorded 
with  the  observation  of  each  failure,  such  that,  comparisons  could  be  made  with  the 
expected  failures  generated  by  the  failure  rate  models. 
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4.  Summary  Statistics 

Table  6  summarizes  the  two  data  sets  and  provides  such  information  as  the 
number  of  aircraft  assigned,  sorties  flown,  flight  hour  totals,  component  failures,  etc. 


TABLE  6 

SUMMARY  STATISTICS 

DATA  SET 

AIRCRAFT    SORTIES 

FLT  HRS  SORTIE 
FLT  HRS      MEAN      STD.DEV. 

MODEL 

30                   3489- 

7608.8            2.18             .98 

VALIDATION 

44                  3012 

7030.4            2.33             .54 

WUC 

MODEL  SET 
FAILURES 

VALIDATION  SET 
FAILURES 

734H100 

105 

58 

46X1600 

34 

49 

4622100 

49 

6 

56X2500 

85        • 

79 

632Z100 

49 

57 

74A1500 

107 

75 

74A5M00 

118 

104 

5772200  . 

34 

33 

6918100 

81 

93 

MODEL  SET  = 

VF 

-I14and  VF-213 

VALIDATION  SET 

=  VF-11,  VF-14,  VF-31,and  VF-32 
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III.  FLIGHT  HOUR  MODELS 

A.  OVERVIEW 

The  two  models  discussed  in  this  chapter  describe  component  failures  as  a 
function  of  flight  hours  only.  The  first  fiight  hour  model  to  be  examined  is  the 
exponential,  which  is  the  model  currently  used  by  the  Navy  to  describe  the  distribution 
of  the  time  between  component  failures.  This  model  has  some  weaknesses  inherited 
from  the  properties  of  the  exponential  distributio.n.  The  effect  these  problem  areas 
might  have  on  the  model's  ability  to  accurately  predict  component  failures  will  be 
discussed  in  the  next  section.  The  three  models  offered  as  alternatives  should  provide 
improvement  to  some  or  all  of  these  problem  areas,  and  the  comparisons  made  in 
Chapter  V  should  shed  new  light  on  the  magnitude  of  these  deficiencies. 

The  second  fiight  hour  model  discussed  is  the  WeibuU.  It  is  offered  as  an 
alternative  method  of  describing  the  time  between  failures  when  fiight  hours  are  the 
only  factors  contributing  to  the  component  malfunction. 

B.  EXPONENTIAL  MODEL 

1.  Weaknesses 

There  are  two  major  weaknesses  associated  with  this  model  that  cause  some 
concern.  First,  the  exponential  is  a  fiight  hour  model  and  as  such  it  is  insensitive  to  the 
effects  of  the  sortie  on  the  component's  time  of  failure.  In  Chapter  V  it  will  be  shown 
that  the  expected  number  of  component  failures  for  the  time  interval  (O.t]  is 
proportional  to  t,  so  that  the  expected  number  of  failures  for  a  specific  fiight  hour  total 
is  not  affected  by  the  average  sortie  length  (i.e.,  the  number  of  sorties  fiown). 

Second,  the  exponential  distributions  are  "NO  WEAR"  distributions,  which 
indicates  that  a  component's  probability  of  surviving  a  fiight  of  duration  t  is  the  same 
for  new  components  as  it  is  for  components  which  are  "S"  hours  old  {i.e..  the 
conditional  survival  probability  is  equal  to  the  unconditional  probability  of  survival). 
This  is  the  "memor\'less"  property  of  the  exponential  distribution,  and  implies  that  the 
prior  use  of  a  component  has  no  effect  on  the  time  of  its  failure. 

The  same  concern  is  refiected  in  the  exponential  distribution's  constant  hazard 
rate  r(t)  (i.e..  r(t}=}v  for  all  values  of  t).  As  such,  the  failure  rate  X  would  not  be 
affected  by  the  time  on  test  t.  It  is  believed  that  certain  components  exhibit  hazard 
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rates  that  are  decreasing  functions  of  time  (t).  In  fact,  many  pieces  of  electronic 
equipment  are  routinely  exposed  to  burn-in  periods  to  make  the  component  less 
susceptible  to  failure  when  put  into  operation.  Conversely,  it  can  be  argued  that  other 
components  have  hazard  rates  that  increase  with  time  (i.e.,  a  component  is  more  likely 
to  fail  as  it  experiences  some  wear  and  tear).  The  exponential  distribution  is  not 
capable  of  modelling  these  diflerences. 
2.  Ma.ximum  Likelihood  Results 

The  maximum  likelihood  estimate,  X',  of  the  exponential  failure  rate 
parameter  "k  (failures/ hour)  is  found  by  iteratively  solving  equation  3.1r 

^         _  t.exp(-Xt.)  T:,t:>0,  X>0,  m<n 

•   yT.  =  V_L_J2^ \L.  I'J       '  (eqn3.1) 

Y    •     Y  ^-^^P^-^V  i=U2,...,n    j=l,2 m 

where,  the  left  hand  side  of  equation  3.1  equals  the  aggregate  flight  time  for  those 
sorties  that  did  not  experience  a  specific  component  failure,  t:  is  the  duration  of  the 
sortie  causing  the  failure,  n  is  the  total  number  of  records,  censored  and  uncensored, 
within  the  final  data  set,-  and  m  is  the  number  of  uncensored  observations  (failures). 

If  the  sortie  lengths,  t:.  are  equal  to  t  for  all  j,  then  equation  3.1  can  be 
simplified  to  the  closed  form  expression  given  in  equation  3.2.  Thus,  the  average  sortie 
length  can  be  used  to  provide  a  "rough"  estimate  of  the  maximum  likelihood  estimator 
I': 

-ln(FTWFTFT) 
X'  = (eqn3.2) 


where,  FTWF  (Flight  Time  Without  Failure)  equals  the  sum  of  the  T-'s  for  i=  l,...n, 
TFT  equals  the  total  flight  time  for  the  deployment,  and  t  equals  the  average  sortie 
length. 

The  maximum  likelihood  estimate  for  X'  was  found  using  the  FORTR.A.X 
program  listed  in  Appendix  E.  The  parameter  estimates  and  95  percent  confidence 
limits  for  these  estimates  are  displayed  in  Table  7. 


"For  the  derivation  of  this  maximum  likelihood  equation,  the  reader  is  referred  to 
Appendix  D. 
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TABLE  7 
EXPONENTIAL  MAXIMUM  LIKELIHOOD  RESULTS 

WUC  X'  di'  CL^ 


734H100 

.0140 

.0018 

(.0106,.0175) 

46X1600 

.0048 

.0010 

{.0029,.0067) 

4622100 

.0065 

.0011 

(.0043,.0086) 

56X2500 

.0114 

.0021 

(.0073,.0155) 

632Z100 

.0067 

.0012 

(.0043,.0090) 

74A1500 

.0143 

.0019 

(.0105,.0181) 

74A5M00 

.0159 

.0024  ■ 

(.0113,.0206) 

5772200 

.0048 

.0011 

(.0026,.0070) 

6918100 

.0108 

•  .0015 

(.0079,.0137) 

** 


"*   CL=  95  PERCENT  CONFIDENCE  LIMITS 


3.  Quantile-Quantile  Plots 

A  graphical  method  often  used  to  determine  an  empirical  data's  goodness  of 
fit  with  respect  to  a  theoretical  distribution  is  the  quantile-quantile  plot.  The  technique 
requires  the  observations,  X^,X2,...,X^,  to  be  ordered  from  smallest  to  largest. 
Empirical  quantiles  are  then  computed  for  these  order  statistics  based  on  a  formula, 
such  as  Qj  =  (i-.5);n,  where  Q.  is  the  empirical  quantile  for  the  X/|^  order  statistic  and  n 
equals  the  total  number  of  observations.  These  empirical  quantiles  are  plotted  against 
the  theoretical  quantiles  to  produce  a  graphical  estimate  of  goodness  of  fit. 

Some  adjustments  to  the  above  technique  must  be  made  to  overcome  the 
problems  of  censoring.  The  first  adjustment  relates  to  the  problem  associated  with 
incomplete  knowledge  of  the  actual  time  of  component  failure.  In  Chapter  II,  it  was 
mentioned  that  the  failure  occurred  during  the  interval  T.  to  T. +  t..  To  simplify  this 
problem,  it  is  assumed  that  the  midpoint  of  that  interval  represents  the  actual  time  of 
failure. 
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The  second  adjustment  that  needs  to  be  made  accounts  for  those  observations 
that  are  right  hand  censored.  The  Kaplan-Meier  method  for  quantile  estimation  will 
account  for  those  observations  and  will  provide  a  more  accurate  estimate  for  the 
empirical  quantiles.  This  method  also  requires  that  each  observation,  censored  and 
uncensored,  be  ordered  from  smallest  to  largest.  The  empirical  quantiles  (Q.)  are 
computed  using  equation  3.3: 

Qj=l-[(n  +  .5)/n]  n  [(n-k+.5)/(n-k+1.5)]      for  i  in  I  (eqn  3.3) 

kin  I 
k<i 

where,  I  is  the  set  of  all  uncensored  observations,  and  n  is  the  total  number  of 
observations,  censored  and  uncensored.   [Ref  7:  p.  234] 

An  exponential  quantile-quantile  plot  was  constructed  for  each  of  the  nine 
components  studied.  The  quantile-quantile  plots  for  work  unit  codes  734H100  and 
56X2500  are  displayed  in  Figures  3.1  and  3.2.  and  the  remaining  quantile-quantile  plots 
are  presented  in  Appendix  F.  The  y-axis  of  the  plot  represents  the  quantiles,  P(T:^t), 
and  the  x-axis  represents  the  time  between  component  failures.  The  asterisks  identify 
the  uncensored  observations,  the  x  's  identify-  the  locations  of  the  censored 
observations,  and  the  theoretical  distribution  evaluated  using  X'  is  defined  by  the  solid 
line.  The  empirical  data  is  fit  well  by  the  theoretical  distribution  if  the  data  points  are 
tightly  packed  around  the  solid  hne.  Data  points  which  deviate  from  the  line  can 
provide  information  about  areas  of  weakness. 

Figure  3.1  illustrates  the  best  exponential  fit  of  all  the  quantile-quantile  plots, 
but  the  plot  still  shows  the  lower  to  mid  quantiles  associated  with  the  empirical  data 
points  overestimate  the  quantiles  computed  using  the  theoretical  distribution.  Figure 
3.2  provides  a  better  representation  of  the  other  quantile-quantile  plots.  In  even.'  plot, 
the  empirical  quantiles  overestimate  the  theoretical  quantiles  for  quantiles  less  than  .5, 
and  underestimate  the  theoretical  quantile  for  quantiles  greater  than  .75. 

Figure  3.3  displays  the  WeibuU  cumulative  distribution  function  for  a  equal  to 
.5,  1,  and  2.  The  reader  should  be  aware  that  a  WeibuU  distribution  with  shape 
parameter  (a)  equal  to  1  is  the  exponential  distribution.  Figure  3.3  shows  that  for 
values  of  time  less  than  tg,  the  cumulative  distribution  for  the  curve  represented  by 
a  =.5  is  greater  than  the  exponential  curve  (a=  1),  and  for  values  of  time  greater  than 
tg  it  is  less   than  the  cumulative  distribution  of  the  exponential  curve.  This  is  the 
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Figure  3.1     Exponential  Quantile-Quanlile  Plot:  WUC  734H100. 
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Figure  3.2     F:xponcntial  Quantile-Quantile  Plot:  WL'C  56X2500. 
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identical  situation  displayed  in  each  of  the  nine  quantile-quantile  plots.  This  would 
suggest  that  the  distribution  representing  the  time  between  component  failures  might 
be  better  modelled  with  a  WeibuU  distribution  with  parameters,  X>0  and  0  <  a  <  1. 


Figure  3.3     Weibull  Cumulative  Distribution  Function. 

C.       WEIBULL  MODEL 

1.  Properties 

The  next  flight  hour  model  that  will  be  discussed  is  the  Weibull  distribution, 
which  can  be  a  viable  alternative  to  the  exponential  distribution  since  it  is  capable  of 
modelling  many  different  failure  rates.  The  major  diflerence  between  the  Weibull  and 
Exponential  models  is  that  the  former  is  usually  not  "memor\'less",  which  indicates  that 
the  component's  previous  usage  can  be  exploited  in  the  determination  of  the  expected 
failure  time. 
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The  two  parameters  associated  uith  the  WeibuU  distribution  are  the  scale 
parameter  (X)  and  the  shape  parameter  (a),  with  each  shape  parameter  defining  a 
specific  failure  rate  distribution.  As  such,  the  WeibuU  family  of  distributions  has  the 
capability  to  model  many  difierent  types  of  component  wear.  For  example,  a  values 
greater  than  one  describe  a  family  of  "WEAR  OUT"  distributions  (increasing  failure 
rate  functions);  an  a  value  equal  to  one  describe  the  "NO  WEAR"  exponential:  and  a 
values  less  than  one  describe  a  family  of  "WEAR  IN"  distributions  (decreasing  failure 
rate  functions).  Consequently,  the  WeibuU  distribution  provides  an  additional 
modeUing  fiexibiUty  that  should  reduce  the  latter  two  areas  of  concern  associated  with 
the  exponential  distribution. 

The  WeibuU  cumulative  distribution  function  (Figure  3.3)  and  the  exponential 
quantile-quantile  plots  indicate  that  a  WeibuU  model  with  scale  parameter  (X>0)  and 
shape  parameter  (a<l)  might  provide  appropriate  descriptors  of  the  failure  rate 
distribution  for  each  of  the  components  studied.  As  mentioned  above,  a  shape 
parameter  less  than  one  is  indicative  of  a  "WEAR  IN"  family  of  distributions. 
Distributions  of  this  type  are  characterized  b\  a  hazard  rate  function.  r(t).  which  is 
decreasing  in  t,  (i.e.,  the  component's  failure  rate  decreases  as  t  increases).  Since,  the 
majority  of  components  studied  are  electronic  components  (AVIONICS),  this 
assumption  of  a  decreasing  failure  rate  is  probably  an  improvement  over  the  constant 
failure  rate  of  the  exponential  model. 
2.  Maximum  Likelihood  Results 

Maximum  Ukelihood  solutions  can  be  obtained  by  setting  the  partials  equal  to 
zero,  and  solving  for  X'  and  a'  using  nonlinear  optimization  techniques  similar  to  the 
one  described  in  Appendix  G.  The  partial  derivatives  of  the  log-likelihood  function 
with  respect  to  the  parameters  X  and  a  are  displayed  belowr 

BK_^      aX«-'|(Ti4-ti)«-Ti«|ST-(ti)     y^^a-l-j-a 

dX       \  l-S-pjt;)  -  ) 


The  derivation  of  these  equations  is  given  in  Appendix  D. 
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where,  a,X>0;  T-,t^,Tj^O;  i=  1,2 m  (the  set  of  uncensored  observations-  failures); 

and  j  =  m+  l,m+2,...,n  (the  set  of  censored  observations). 

The  maximum  likelihood  estimates  for  X'  and  a'  were  found  using  the 
FORTRAN  program  listed  in  Appendix  G.  The  parameter  estimates  and  95  percent 
confidence  limits  for  these  estimates  are  displayed  in  Table  8. 


TABLE  8 

WEIBULL  MAXIMUM  LIKELIHOOD  RESULTS 

wuc 

X' 

CLx 

a' 

CL„ 

734H100 

.0149 

(.0084,.0216) 

"^  .8079 

(.8014,.8146) 

46X1600 

.0044 

{.0006,.0082.) 

.6664 

(.6627..6703) 

4622100 

.0066 

(.0024,.0108) 

.7708 

(.7666,.7750) 

56X2500 

.0119 

(.0063,.0176) 

.8170 

{.8114,.8227) 

632Z100 

.0067 

(.0026,.0109) 

.7970 

(.7928,.8011) 

74A1500 

.0165 

•     (.0084,.0246) 

.7066 

{.6985,.7147) 

74A5M00 

.0164 

(.0104,.0224) 

.8844 

(.8784,.8904) 

5772200 

.0045 

(.0009,.0082) 

.7385 

(.7349,.7422) 

6918100 

.0112 

(.0058,.0167)' 

.8142 

(.8087,:8196) 

***   CL  = 

95  PERCENT  CONFIDENCE  LIMITS 

3.  Quantile-Quantile  Plots 

The  parameter  estimates  listed  in  Table  8  were  used  to  construct  Kaplan- 
Meier  quantile-quantile  plots  for  each  of  the  nine  components  studied.  Two  of  these 
quantile-quantile  plots  are  displayed  in  Figures  3.4  and  3.5,  and  the  other  seven  are 
included  in  Appendix  H.  Both  Figures  graphically  illustrate  that  the  data  is  fit  well  by 
the  Weibull  distribution,  with  the  worst  case  fit  represented  by  Figure  3.4.  This  plot 
shows  that  the  Weibull  distribution  has  a  tendency  to  underestimate  the  empirical 
quantiles  below  ten  percent,  but  fits  the  data  extremely  well  for  quantiles  greater  than 
ten  percent. 
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The  results  of  the  Weibull  quantile-quantile  plots  are  not  conclusive,  but  offer 
credibility  to  the  use  of  the  Weibull  model  as  an  alternative  method  of  modelling 
component  failure  rates  when  flight  time  is  the  only  contributing  factor. 
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Figure  3.4     Weibull  Quantile-Qiiantile  Plot:  WUC  734H100. 
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Figure  3.5     Weibull  Quantilc-Quantile  Plot:  WLC  56X2500. 
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IV.  SORTIE  DEPENDENT  MODELS 

A.  OVERVIEW 

The  two  models  that  will  be  discussed  in  this  chapter  differ  from  the  models  in 
Chapter  III  in  that  component  failures  can  be  affected  by  stresses  incurred  as  a  result 
of  the  sortie.'^  These  stresses  refer  to  events  that  occur  during  each  flight,  regardless  of 
the  sortie  length,  and  would  include  events  such  as  equipment  initialization,  takeoffs, 
landings,  etc.  The  first  sortie-dependent  model  to  be  discussed  is  the  Geometric  model, 
which  determines  a  constant  probability  that  a  component  fails  on  a  given  sortie. 

The  second  sortie-dependent  model  that  will  be  briefly  looked  at  is  the  Vlixed 
model  described  in  reference  5.  This  model  could  also  be  classified  as  a  flight  hour 
model  since  it  determines  the  probability  of  a  component  failure  based  on  both  sortie 
and  flight  hour  contributing  factors.  It  is  incorporated  in  this  section  because  it  also 
utilizes  the  geometric  distribution  to  determine  the  probability  that  the  first  failure 
occurs  on  the  s  sortie.  The  major  difference  between  these  two  models  is  that  the 
Mixed  model's  probability  of  failure  for  a  given  sortie  is  not  constant.  In  fact,  it  is  a 
function  of  a  constant  failure  probability  caused  by  sortie-related  stress,  and  a 
probability  of  failure  which  is  dependent  on  the  sortie  length. 

B.  GEOMETRIC  MODEL 

1.  Properties 

The  basic  premise  of  the  Geometric  model  is  that  the  component  failure 
results  due  to  sortie-related  stress  and  is  insensitive  to  differences  in  sortie  length.  This 
assumption  can  then  be  used  to  formulate  the  component's  survival  or  failure  for  a 
given  sortie  as  a  Bernoulli  trial,  where,  the  component  fails  with  probability  p  or 
survives  with  probability  1-p.  The  geometric  distribution  is  used  to  determine  the 
probability  of  observing  the  first  failure  on  the  s^^  Bernoulli  trial. 

The  geometric  distribution  is  in  many  ways  similar  to  the  exponential 
distribution  and  might  be  considered  the  exponential's  discrete  analog.  Like  its 
continuous  counterpart,  the  geometric  distribution  is  "memor\'less"  and  has  a  constant 
failure  rate,  which  would  indicate  that  the  Geometric  model  also  ignores  the  history  of 


A  sortie  refers  to  one  complete  flight  evolution,  take-off  to  full-stop  landing. 
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component  usage.  Any  improvement  ofTered  by  this  model  would  be  attributed  to  the 
efiects  of  the  sortie  on  the  component's  failure. 

It  should  be  apparent  that  a  low  variability  in  sortie  length  (i.e.,  the  sortie 
length  is  approximately  constant)  would  imply  that  the  exponential  distribution  is 
equivalent  to  a  geometric  distribution  with  p=  l-exp(-)tt);  where  t  equals  the  average 
sortie  length.  In  this  case,  both  models  would  provide  similar  estimates  of  component 
failures.  The  flight  data  used  to  generate  the  parameter  estimates  was  acquired  from 
aircraft  flying  during  Flex-deck  operations,  and  should  have  greater  variability  than 
would  be  expected  from  aircraft  flying  under  normal  cyclic  operations.^ 
2.  Maximum  Likelihood  Results 

The  maximum  likelihood  expression  for  p',  the  estimate  of  the  constant  sortie 
failure  probability  is  given  in  equation  4.1:^ 

m 
^'  =  8:=  1,2,...     i=1.2 n  (eqn4.1) 


y  s 


where,  m  represents  the  number  of  uncensored  observations  (sorties  with  component 
failures),  and  the  summation  of  the  s-'s  for  i=  1.2,...,n  represents  the  total  number  of 
sorties  flown.  Thus,  the  maximum  likelihood  estimator  p'  equals  the  number  of  flights 
experiencing  a  component  failure  divided  by  the  total  number  of  sorties  flown.  This  is 
the  same  maximum  likelihood  estimate  which  would  have  resulted  without  censoring. 
The  maximum  likelihood  estimates  p'  for  p  and  the  associated  95  percent  confidence 
intervals  are  listed  in  Table  9. 

C.       MIXED  MODEL 

1.  Properties 

This  section  will  provide  a  brief  description  of  the  Mixed  model  proposed  by 
Le\7  [Ref  5].  This  particular  model  is  structured  such  that  the  sortie  and  the  sortie 
length  (flight  hours)  can  combine  to  cause  a  component  failure.  It  is  obvious  that 
factors  contributing  to  one  component's  failure  may  or  may  not  be  detrunental  to 
another  component's  operating  life.  Most  likely,  "the  life  of  electronic  components  is 
probably  affected  more  by  the  number  of  power  surges  to  which  they  are  subjected 


"^Duruis  Flex-deck  operations,  aircraft  sortie  length  is  not  constrained  to  a  one  to 
two  hour  launch  and  recover."  cycle. 

^The  derivation  of  this  equation  is  given  in  Appendix  D. 
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than  operating  hours;  landing  gear  life  should  be  influenced  primarily  by  takeoffs  and 
landings;  and  engine  life  should  be  affected  by  rotational  speed  changes,  which  are 
accompanied  by  heat  transients  and  pressure  changes"  [Ref  8:  p.  2].  The  Mixed  model 
appears  to  have  an  increased  flexibilty  to  model  some  of  these  sortie  and  flight  hour 
factors. 

As  in  the  previous  section,  a  constant  probability  of  failure  is  employed  to 
model  the  effects  of  sortie  induced  stress.  In  addition,  the  Mixed  model  incorporates  a 
failure  probability  for  those  components  whose  lifetime  is  affected  by  the  flight  hours 
fiowTi.  Combining  the  constant  sortie  failure  probability  with  the  flight  hour  probability 
results  in  the  expression  for  the  probability  of  failure  p-  for  the  i  sortie,  which  is 
expressed  as  the  complement  of  a  component  surviving  the  stresses  of  the  sortie  and 
the  flight  hours  flown: 

Pi=  l.exp(-^)(l.pQ)       0<pQ<l,X>0,t->0  (eqn4.2) 

where,  the  flight  hour  failure  distribution  is  exponential  and  p  equals  the  constant 
sortie  failure  probability.    - 

2.  Maximum  Likelihood  Results 

The  partial  derivative  of  the  log-likelihood  function  with  respect  to  the 
parameters  X  and  p^  are  listed  below:^ 

^K      ^     {l-p_)t:exp(-}lt:)  ^ 

—  =  y ^    '      ^  -  T;  -  y  T:  (eqn  4.3) 

dX     X   l-exp(->.t^)(l-pQ)       1    j-    J  ^  4         ; 

^K  exp(->.t:)  S:-I  S: 

—  =  y -— ^ i X ^-  (eqn  4.4} 

^p^     Yl-expl-Xtpd-p^)        1-p^       Y     1-p^ 

where,  0<Pq<1;  X>0;  T-,t-,T:>0;  sj=l,2,...;  i=  l,2,...,m  (the  set  of  uncensored 
observations  -  failures);  and  j  =  m+  l,m+2,....n  (the  set  of  censored  observations). 

Maximum  likelihood  solutions  for  X'  and  p  '  were  obtained  using  the 
FORTR.A.N  program  listed  in  Appendix  G.  For  each  of  the  nine  components  studied. 
the  log-likelihood  function  was  maximized  when  X'  equaled  zero,  and  p  -  equaled  p', 
the  maximum  likelihood  estimate  for  the  constant  failure  rate  p  of  the  Geometric 


The  derivation  of  these  equations  is  given  in  Appendix  D. 
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model.  The  results  for  X'  are  included  in  Table  9.  The  reader  should  note  that  the 
values  for  p'  and  Pq'  are  identical,  and  the  failure  for  X  is  zero  for  every  component. 
Since  both  models  provide  the  same  maximum  likelihood  results,  they  will  be  jointly 
referred  to  as  the  Geometric  model  when  compared  with  the  other  two  models. 


TABLE  9 

MAXIMUM  LIKELIHOOD  RESULTS 
SORTIE  DEPENDENT  .MODELS 

wuc 

p'** 

CLp 

X'* 

CLX     ■ 

734H100 

.0301 

(.0188,.0414) 

.0000 

(.0000,.0055) 

46X1600 

.0097 

{.0038, .0156) 

.0000 

(.0000,.0027) 

4622100 

.0140 

{.0069,.0211) 

.0000 

(.0000,.0034) 

56X2500 

.0244 

(.0147,.0341) 

.0000 

{.0000,.0047) 

'      632Z100 

.0140 

(.0069,.02ll) 

.0000' 

{.0000..0033) 

74AI500 

.0307 

{.0191,.0423) 

.ooco 

(.0000,.0056) 

74A5M00 

.0338 

(.0228,.0448) 

.0000 

(.0000,.0053) 

5772200 

.0097 

(.0037,.0l57) 

.0000 

(.0000,.0028) 

■    6918100 

.0232 

(.0137,.0327) 

.0000 

{.0000,.0045) 

CL=  95  PERCENT  CONFIDENCE  LIMITS 

*      MIXED  MODEL  ONLY 

**    Po'  =  P' 
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V.  MODEL  COMPARISONS 

A.  PURPOSE 

The  purpose  of  this  chapter  is  to  compare  the  Exponential  model  with  the 
Weibull  and  Geometric  models  to  determine  if  the  distribution  currently  used  provides 
the  most  accurate  description  of  the  failure  process.  As  a  method  of  comparison,  the 
Chi-square  goodness-of-fit  test  was  not  considered  for  two  reasons.  First,  the  failure 
rate  information  provided  by  the  censored  observations  would  have  been  lost  because 
of  the  impossibihty  of  determining  the  exact  time  of  component  failure.  Second,  it  was 
considered  more  important  to  examine*  each  model's  capability  to  predict  the  number 
of  component  failures  resulting  from  a  specified  number  of  sorties  and  flight  hours 
flown.  The  next  section  describes  the  methods  of  estimating  failures  for  each  model. 
These  estimates  are  then  compared  with  the  actual  number  of  failures  to  determine 
each  model's  goodness-of-fit. 

B.  COMPONENT  FAILURE  ESTIMATION 

1.  Geometric  Model 

The  determination  of  the  expected  number  of  component  failures  -associated 
with  the  Geometric  model  is  straightforward.  The  parameter  p  estimated  using  the 
Geometric  model  defines  the  probability  that  a  component  would  fail  during  a  given 
sortie.  Furthermore,  if  it  is  assumed  that  the  outcome  of  a  sortie  is  independent  of  the 
outcome  of  the  other  N-1  sorties,  and  p  is  constant  for  each  sortie,  then  the  probability 
law  defining  the  number  of  component  failures  for  the  deployment  is  the  sum  of  N 
independent  Bernoulli  trials.  This  random  sum  has  a  Binomial  distribution  with  an 
expected  value  equal  to  N  x  p.  Thus,  using  the  Geometric  model  the  expected  number 
of  component  failures  is  calculated  as  the  product  of  the  number  of  sorties  flown  (N) 
and  the  probability  that  the  component  fails  during  the  sortie. 

2.  Exponential  Model 

To  compute  the  expected  number  of  component  failures  using  the  Exponential 
model,  each  failure  can  be  thought  of  as  an  occurrence  within  a  Renewal  process;  a 
stochastic  process  which  counts  the  number  of  occurrences  within  an  interval  {0,t].  It  is 
assumed  that  the  times  between  successive  occurrences  are  independently  and 
identically  distributed  non-negative  random  variables.  As  such,  it  can  be  shown  that 
the  expected  number  of  occurrences  within  an  interval  (0,t]  is  equal  to: 
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'^(0  =  lFi,(t)  k=  1,2,3,... 

k 

where,  M(t)  equals  the  expected  number  of  occurrences  (failures)  within  the  interval 
(0,t],  and  Fj^(t)  is  the  k-fold  convolution  of  the  cumulative  distribution  function.  The 
derivative  of  M(t)  with  respect  to  t  defines  the  renewal  rate  of  the  process  as: 

M'(t)  =  m(t)  =  Xrk(^)       ^=12,3,... 
k 

Since  the  k-fold  convolution  of  the  exponential  probability  density  function  is  known 
to  be  Gamma,  m(t)  is  equal  to  the  infinite  sum  of  Gamma  probability  density  functions 
with  scale  parameter  X  and  shape  parameter  k. 

^k^k-l,-Lx                       (^^^k-l 
m(t)  =  V =  Xq-^^  y k=  1,2,3,... 

t       r(k)  J-  r(k) 

In  the  equation  above,  Y  [(X,x)      ,T{k)]  for  k=  1.2.3,...  is  easilv  recognized  as  the 
Taylors  senes  expans,on"of  e^^.    Thus,  the  renewal  rate  for  the  counting  process 
.defined  by  exponential  interoccurrence  times  is  X.    The  expected  number  of" failures 
M(t)  is  found  by  integrating  m{t)  from  0  to  t. 

t 

M(t)  =  JX  ds  =  >.t 
0 

Therefore,  the  expected  number  of  component  failures  is  equal  to  the  product  of  the 
failure  rate  X  and  the  time  on  test  t.  This  particular  Renewal  process  is  referred  to  as  a 
Poisson  process. 

3.  Weibull  Model 

Estimating  the  number  of  expected  failures  associated  with  the  Weibull  model 
is  not  as  easy  as  the  other  methods  discussed  above.  The  observation  of  component 
failures  is  still  a  Renewal  process,  but  the  interoccurrence  times  are  distributed 
according  to  a  more  complex  distribution  function.  Finding  the  k-fold  convolution  of 
either  the  cumulative  distribution  function  or  the  probability  density  function  is  an 
extremely    dilTicult    if  not    impossible    task.     For    processes    defined    by    continuous 
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interevent  times  with  finite  means  and  variances,  the  Elementary  Renewal  Theorem 
states  that  M{t)  is  approximately  Normal  with  mean  =  t/n  and  variance  =t(y^;H-^  as  t 
approaches  infinity  [Ref,  9:  p.  290].  Unfortunately,  this  theorem  does  not  help  estimate 
failures  when  t  is  small.  As  an  alternate  method  to  construct  M(t)  when  t  is  small,  a 
simple  computer  simulation  is  used. 

The  inputs  required  for  the  simulation  included  the  number  of  aircraft 
assigned  to  the  carrier  air  wings,  the  Weibull  shape  and  scale  parameters,  and  an 
approximate  flight  hour  distribution.  This  information  was  used  to  repeatedly  simulate 
the  air  wing  operations  during  the  deployment. 

The  aircraft  were  chosen  at  random  to  fly  a  fictitious  mission  of  some  sortie 
length  t.  For  the  modelling  set,  there  were  30  aircraft  assigned.  The  model  assumed 
that  each  aircraft  had  an  equally  likely  (1/30  percent)  probability  of  being  chosen  to  fly 
a  mission.  Realistically,  aircraft  are  not  chosen  at  random  to  fly  missions,  and  most 
hkely  those  aircraft  recently  flown  are  better  prepared  for  future  flights.  A  weighting 
scheme  could  have  been  used  to  increase  the  probabilities  of  those  aircraft  flown 
recently,  and  decrease  the  probability  for  those  aircraft  that  were  undergoing  simulated 
maintenance  activities. 

The  duration  of  the  mission  was  determined  by  a  random  variable  generated 
from  a  Normal  distribution  with  the  mean  and  variance  estimated  by  the  Model  data 
set  (see  Table  6).  Figure  5.1  depicts  a  histogram  of  the  sortie  lengths  for  the  Model  set 
with  an  overlay  of  the  appropriate  Normal  distribution.  Between  the  values  of  Q  to  6 
hours,  the  Normal  distribution  fits  the  histogram  adequately  enough  to  justify  its  use 
as  the  sortie  length's  parent  distribution.  Values  of  randomly-generated  sortie  lengths 
less  than  .5  hours  or  greater  than  6  hours  were  truncated  to  .5  hours  and  6  hours, 
respectively. 

After  the  selection  of  aircraft,  the  conditional  probability  of  failure  was 
computed  for  a  flight  of  sortie  length  t.  The  only  information  required  for  this 
computation  was  the  component's  Weibull  shape  and  scale  parameters,  and  the 
aircraft's  flight  histor\'  since  the  last  component  failure.  The  conditional  probability  of 
failure  could  then  be  calculated  as  described  in  Appendix  D.  The  value  for  the 
conditional  probability  of  failure  was  then  compared  against  a  Uniform  (0,1)  random 
variable,  and  if  the  uniformly-generated  outcome  was  less  than  or  equal  to  the 
conditional  probability  of  failure,  a  component  failure  was  simulated. 


45 


4  6 

SORTIE  LENGTH 


Figure  5.1     Sortie  Lengtli  Histogram  vs  Normal  Distribution. 

The  number  of  component  failures  was  recorded- at  500  hour  intervals  until 
the  completion  of  the  deployment.  The  deployment  was  simulated  100  times,  and  the 
average  number  of  component  failures  was  computed  for  each  500  hour  interval.  These 
averages  were  the  approximations  to  the  expected  number  of  component  failures  M(t) 
at  500,1000,1500,...  hours  etc. 

C.       COMPARISONS 

1.  Methodology 

The  predicted  failures  were  compared  with  the  actual  failures  as  they  occurred 
in  each  of  the  two  deployments  (IJSS  KENNEDY  and  LSS  INDEPENDENCE).  At 
this  point,  it  is  important  to  clarify  the  specific  calculations  of  the  expected  failures.  As 
described  earlier,  the  accumulated  fiight  time  and  total  sorties  flown  were  recorded  with 
the  observation  of  each  failure.  Ihcrefore.  the  exponential  and  geometric  expected 
lailurcs  could  be  computed  Ubing  the  formulas  listed  in  this  chapter.  Unfortunately,  the 
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expected  number  of  WeibuU  component  failures  was  not  always  defined  at  the 
accumulated  flight  time  of  the  actual  failure.  To  provide  such  an  estimate,  it  was 
necessary  to  interpolate^  between  the  WeibuU  expected  failures  associated  with  the 
simulated  failure  times  above  and  below  the  actual  time  of  component  failure.  The 
difference  between  the  expected  and  actual  failures  was  computed  at  the  time  of  each 
component  failure.  These  differences  were  then  summed  and  divided  by  the  total 
number  of  component  failures  to  provide  an  estimate  of  the  average  prediction  error. 
An  average  difference  close  to  zero,  with  small  standard  deviation,  indicates  that  the 
model  is  a  good  predictor  of  actual  failures.  A  negative  average  shows  that  the  model 
has  a  tendency  to  underestimate  the  number  of  component  failures,  and  a  positive 
average  is  indicative  of  overestimation.  Tables  10  and  11  provide  the  average 
differences  and  corresponding  standard  deviations  for  the  "model"  and  "validation"  data 
sets  respectively. 

The  model  providing  the  best  fit  for  a  specific  component  was  chosen  based 
on  the  average  difierence  between  the  predicted  failures  and  the  actual  failures.  Table 
12  provides  a  summary  listing- of  the  optimal  distributions  for  each  component  in  both 
data  sets.  The  reader  is  cautioned  that  the  optimal  was  chosen  based  on  an  average 
difference  closest  to  zero,  and  in  some  cases  the  expected  number  of  failures  was 
adequately  described  by  other  distributions  as  well  as  the  one  listed  as  optimal.  If  two 
or  more  models  provided  similar  estimates,  the  one  with  the  smallest  standard  deviation 
was  chosen  as  optimal.  The  actual  failures  were  plotted  against  the  predicted  failures  to 
further  illustrate  a  distribution's  ability  or  disability  to  estimate  failures.  These  graphs 
are  included  in  Appendix  I. 
2.  Results 

Before  proceeding  with  the  comparisons,  a  caveat  should  be  placed  on  the 
results.  The  ability  of  the  models  to  predict  the  failures  observed  during  the  deployment 
of  the  "validation"  data  set  is  contingent  on  the  accuracy  of  the  model  parameter 
estimates.  For  a  valid  comparison,  it  is  necessary  that  both  deployments  come  from  the 
same  population.  That  is  to  say,  the  flight  operations  were  conducted  under  similar 
conditions.  The  reader  should  be  aware  that  the  "model"  set  aircraft  and  the 
"validation"  set  aircraft  operated  during  different  seasons  and  in  different  geographic 
locations.  The  parameters  estimated  from  the  "model"  set  may  or  may  not  provide  the 


A  quadratic  interpolation  of  the  form  aFT^  +  bFT  +  c=  EF  was  used.  Where.  FT 
equaled  accumulated  fliizht  time  of  failure.  EF  equaled  expected  failures,  and  a.b.c  were 
parameters  estimated  using  the  known  WeibuU  times  and  expected  failures. 
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TABLE  10 

MODEL  SET  AVERAGE  DIFFERENCES 

WUC:734H100 
Average: 
St.Dev. : 

Exponential 

-0.31 
2.58 

Geometric 

1.51 
3.53 

Weibull 

3.20 
2.89 

WUC:46X1600 
Average : 
St.Dev. : 

0.94 
1.91 

0.50 
1.27 

3.65 
1.93 

WUC:4622100 
Average: 
St.Dev. : 

-2.87 
2.54 

-1.91 
1.97 

0.56 
2.52 

WUC: 56X2500 

Average: 
St.Dev. : 

-1.37 
1.94 

0.17 

1.18 

1.43 
1.76 

WUC:632Z100 
Average: 
St.Dev. : 

-0.04 
1.85 

0.08 
2.12 

2..  74 
2.29 

WUC:74A1500 
Average : 
St.Dev.: 

-1.75 
4.77 

-0.03 
3.92 

5.78 
5.69 

WUC:74A5M00 
Average : 
St.Dev. : 

•       -0.53 
3.31 

1.05 
4.45 

0.09 

4.14 

WUC:5772200 

Average : 
St.Dev. : 

-1.49 
2.13 

-1.84 
1.29 

1.95 
2.44 

WUC:6918100 
Average : 
St.Dev. : 

-4.36 
3.23 

-2.97 
2.53 

-1.75 
3.03 

true  estimate  of  the  population  parameters  required  to  accurately  predict  failures  within 
the  "validation"  data  set. 

Within  the  "model"  set,  all  of  the  models  predict  failures  adequately.  This  is 
not  unexpected  because  the  parameters  dcfming  the  distribution  were  estimated  using 
the  fliaht  data  from  the  "model"  set.  Table   1 1   shows  that  this  is  not  true  for  the 
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TABLE  11 

VALIDATION  SET  AVERAGE  DIFFERENCES 

WUC:734H100 
Average : 
St.Dev. : 

Exponential 

20.76 
11.01 

Geometric 

16.44 
8.85 

Weibull 

26.22 

12.01 

WUC:46X1600 

Average : 
St.Dev. : 

-5.31 
4.04 

-8.04 
5.32 

-0.37 
4.17 

WUC:4622100 
Average: 
St.Dev. : 

16.20 
15.33 

14.57 
13.95 

20.90 
16.86 

WUC: 56X2500 

Average : 
St.Dev. : 

-4.23 
4.01 

-7.41 
5.15 

0.48 
3.49 

WUC:632Z100 
Average: 
St.Dev. : 

-5.88 
3.14     ■ 

-8.43 
4.45 

-1.65 
2.42 

. 

WUC:74A1500 
Average : 
St.Dev. : 

5.82 
6.10 

2.00 
4.26 

17.12 
8.06 

WUC:74A5M00 
Average : 
St.Dev. : 

8.67 
3.41 

2.83 
3.75 

10.39 
3.74 

WUC: 5772200 
Average : 
St.Dev. : 

-4.61 
1.93 

-6.33 
2.40 

-0.07 
2.02 

WUC:6918100 
Average: 
St.Dev. : 

-2.01 
7.12 

-5.81 
8.48 

2.57 
7.27 

deployment  defined  by  the  "validation"  data  set.  In  most  cases,  a  single  optimal  model 
can  be  identified. 

Table  10  illustrates  an  extremely  minor  diflerence  between  the  e.xponential  and 
geometric  estimates.  In  fact,  it  could  be  argued  the  geometric  and  the  exponential  were 
interchangeable  as  predictors  of  failures  for  the  "model"  set  deployment.     For  the 
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TABLE  12 

OPTIMAL  DISTRIBUTIONS 

wuc 

MODEL  SET 

VALIDATION  SET 

734H100 

Exponential 

Geometric 

46X1600 

Geometric 

WeibuU 

4622100 

WeibuU 

Geometric 

56X2500 

Geometric 

WeibuU 

632Z100 

Exponential 

WeibuU 

74A1500 

Geometric 

Geometric 

74A5M00 

Geometric 

Geometric 

5772200 

Geometric 

WeibuU 

6918100 

WeibuU 

Exponential 

Inertial  Measuring  Unit  (WUC:734H100),  the  optimal  distribution  is  Usted  as 
exponential  for  the  "model"  set  and  geometric  for  the  "validation"  set,  but  it  could  have 
been  geometric  in  both  cases.  During  the  deployment  represented  by  .the  "vaUdation" 
data  set  the  average  sortie  length  increased  by  approximately  S.5  percent.  As 
mentioned  earlier,  the  exponential  distribution  is  insensitive  to  changes  in  the  average 
sortie  length  and  would  provide  an  estimate  based  on  flight  hours  only.  On  the  other 
hand,  if  the  geometric  distribution  was  correct  the  longer  sortie  lengths  would  imply 
that  a  specific  flight  hour  total  was  accompUshed  with  fewer  sorties  flown.  It  would 
provide  a  smaller  estimate  than  would  the  Exponential  model.  For  Work  Unit  Code 
734H100,  the  Exponential  model  was  unable  to  account  for  the  longer  average  sortie 
length  (in  the  "validation"  set)  and  overestimated  the  actual  number  of  failures. 

The  Geometric  model  appears  to  provide  the  best  fit  for  components  that  are 
not  necessarily  operated  continuously.  Examples  of  components  of  this  type  would 
include  the  AWG-Q  Radar  Transmitter  (WUC:74A1500)  and  Detail  Data  Display 
(WUC:74A5M00).  A  possible  exception  to  this  rule  would  be  the  Inertial  Measuring 
Unit  (WUC:734H100).  Despite  its  continuous  use,  the  Inertial  Measuring  Unit's  gyros 
and  accelerometers  are  extremely  sensitive  to  the  movements  and  impacts  associated 
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with  carrier  takeofis  and  landings,  and  it  is  extremely  likely  that  component  failures  are 
caused  by  sortie-related  stress  vice  continuous  flight  operations. 

For  continuous  use  components  within  the  "validation"  data  set,  the  WeibuU 
distribution  provided  the  best  estimate  of  component  failures.  The  exponential 
distribution  was  the  optimal  predictor  for  only  one  component  (WUC:6918100).  and  in 
that  case,  the  WeibuU  distribution  was  nearly  as  accurate.  From  a  readiness 
standpoint,  the  WeibuU  model  might  be  considered  a  better  choice  because  of  its 
tendency  to  overestimate  the  actual  number  of  failures.  Amazingly,  the  exponential 
distribution  predicted  failures  best  for  less  than  10  percent  of  the  components  within 
the  "vaUdation"  data  set.  In  support  of  the  current  model,  it  should  be  noted  that  the 
exponential's  estimate  of  failures  was  more  than  adequate  for  the  "model"  set. 

As  mentioned  previously,  the  WeibuU  distribution  with  shape  parameter  equal 
to  one  is  exponential.  Figure  5.2  demonstrates  the  capabUity  of  the  WeibuU  simulation 
to  provide  estimates  similar  to  those  generated  for  the  Exponential  model.  This  figure 
shows  the  WeibuU  and  exponential  predictions  for  the  DetaU  Data  Display 
(WUC:74A5M00)  with  estimated  WeibuU  shape  parameter  equal  to  .91. 

It  should  not  be  surprising  that  in  some  cases,  none  of  the  models  will  be 
adequate  estimators  of  the  actual  failures.  This  type  of  phenomenon  was  observed  for 
the  Fuel  Tank  Release  Mechanism  (WUC:4622100).  For  two  successive  deployments, 
total  failures  went  from  49  to  6  for  a  similar  number  of  flight  hours  and  sorties  flown. 
Based  on  the  available  data,  it  is  hard  to  determine  if  49  was  too  high  or  6  was  too 
low.  This  type  of  outcome  emphasizes  the  difficulty  in  predicting  the  actual  failures  for 
an  upcoming  deployment  based  on  parameters  estimated  from  one  deployment  rather 
than  the  entire  historical  data  base. 
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Figure  5.2     Flight  Hour  Models:  Model  Set  WUC:74A5M00. 
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VI.  SUMMARY 

A.       ADVANTAGES  AND  DISADVANTAGES 

This  section  provides  a  brief  summary  of  the  primary  advantages  and 
disadvantages  associated  with  each  of  the  four  models  discussed  in  Chapters  III  and 
IV.  The  major  advantage  of  the  Exponential  and  Geometric  models  is  that  the  number 
of  expected  failures  is  easily  calculated  using  closed  form  expressions.  Both  models 
express  failures  as  a  linear  function  of  only  one  factor,  flight  hours  (Exponential)  or 
sorties  flown  (Geometric).  There  is  a  tendency  then  to  concentrate  on  aggregate  flight 
totals  while  ignoring  the  specific  flight  history  of  each  aircraft.  As  a  result,  neither 
model  is  able  to  adjust  its  failure  estimates  to  compensate  for  differences  in  flight 
programs. 

The  Mixed  model,  on  the  other  hand,  attempts  to  account  for  failures  resulting 
from  sortie  stress  and;  or  continuous  use.  Unlike  the  previous  models,  the  expected 
number  of  component  failures  is  a  function  of  a  specific  flight  program.  The 
probability  of  failure  associated  with  the  Vlixed  model  is  a  function  oi"  the  sortie 
duration.  As  such,  this  model  should  be  more  sensitive  to  difTerences  in  flight* 
programs  and  would  provide  a  more  accurate  estimate  of  failures  than  the  Exponential 
and  Geometric  models.  Use  of  the  Exponential  distribution  to  model  continuous-use 
components  has  weakened  the  forecasting"  capabilities  of  the  Vlixed  model.  Lack  of 
consideration  for  flight  hours  flown  since  the  last  component  failure  oversimplifies  the 
problem.  As  stated  previously,  the  small  variability  in  sortie  length  combined  with  the 
"memoryless"  property  of  the  Exponential  distribution  would  imply  that  the  probability 
of  night  hour  induced  failure  is  approximately  constant  for  all  sorties  flown.  Since  the 
stress  related  failure  probability  is  also  constant,  this  model  would  not  necessarily 
predict  any  better  than  the  Geometric  model.  This  type  of  phenomenon  was  observed 
in  this  study. 

The  WeibuU's  strongest  feature  is  its  ability  to  model  failure  rates  for  components 
which  have  increasing,  decreasing,  or  constant  failure  rates.  Since  it  does  not  possess 
the  memor>'less  property  of  the  Exponential/Geometric  models,  the  expected  number  of 
failures  is  dependent  upon  previous  flight  histories.  Like  the  Mixed  model,  this  should 
provide  an  improved  estimate  of  component  failures.  One  disadvantage  is  a  lack  of 
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simple  mathematical  formulas  with  which  to  calculate  the  expected  failures.    However, 
in  today's  world  of  computers  this  should  not  be  a  major  problem. 

B.       CONCLUSIONS  AND  RECOMMENDATIONS 

It  is  extremely  difficult  to  estimate  the  number  of  component  failures  which 
occur  during  a  deployment  based  solely  on  the  data  provided  by  another  deployment. 
Such  an  approach  would  not  account  for  differences  in  flight  and  maintenance 
personnel,  weather  conditions,  geographical  location,  mission  type,  aircraft's 
maintenance  histor>',  etc.  For  example,  the  number  of  actual  failures  observed  for  the 
Fuel  Tank  Release  Mechanism  went  from  forty-nine  to  six  for  a  similar  number  of 
flight  hours  and  sorties  flown.  To  minimize  the  effects  caused  by  such  variations  in 
deployments  it  is  recommended  that  these  models  be  evaluated  using  an  enlarged  data 
base. 

Surprisingly,  in  this  study,  the  model  currently  used  by  the  Nav\'  (Exponential) 
did  not  provide  the  best  estimate  of  component  failures.    The  Geometric  provided  a 
better  fit  for  components  which  were  not   operated  continuously,  and  the  WeibuU. 
performed  better  when  the  components  were  operated  continuously.    In  fact,  overall, 
the  Exponential  was  the  least  effective  model  for  these  nine  components. 

The  results  obtained,  using  the  Weibull  simulation  were  encouraging."  The  model 
produced  satisfactory'  estimates  with  a  minimal  number  of  assumptions,  data  input,  and 
software  coding.  While  the  results  obtained  using  this  model  were  not  overwhelming, 
an  improvement  in  prediction  capability  was  observed. 

Two  possible  directions  for  improving  the  current  method  of  estimating  demand 
(i.e.,  component  failures)  are  the  construction  of  a  universal  model  and  the  use  of  high 
resolution  simulations.  The  first  alternative  suggests  the  use  of  a  model  which  can 
describe  components  with  different  failure  rate  functions,  and  incorporates  a 
methodology  utilizing  soriie-induced  stress  and  flight  hours  flown  to  compute  an 
estimate  of  component  failures.  The  model  currently  used  does  not  have  this 
capability. 

A  slight  modification  to  the  Mixed  model  described  in  Chapter  IV  might  provide 
the  modelling  flexibility  desired.  A  failure  rate  distribution,  such  as  the  Weibull,  could 
be  used  to  describe  continuous-use  failures  instead  of  the  Exponential.  Such  a  change 
would  ensure  that  the  component's  failure  probability  could  be  influenced  by  prior  use. 
That  is  to  say.  the  probability  of  failure  for  a  given  sortie  would  now  be  a  function  of 
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the  sortie  length  and  the  aggregate  flight  time  since  the  last  component  failure.  The 
strengths  of  the  Weibull  distribution  combined  with  the  constant  probability  of  failure 
(sortie)  associated  with  the  Mixed  model  should  provide  an  estimate  of  demand  which 
is  more  sensitive  to  variations  within  and  between  flight  programs. 

The  price  paid  for  this  improvement  in  forecasting  ability  is  an  increase  in  the 
complexity  of  the  expected  demand  computations.  With  the  abundance  of  computers 
within  the  fleet  this  should  not  be  a  problem.  Perhaps,  ten  to  twenty  years  ago  it  was 
necessary  to  rely  on  simple  models  for  prediction  and  the  crude  estimates  they 
provided.  This  is  not  the  case  today;  computer  simulations  can  be  used  when  simple 
mathematical  formulas  are  not  available.  As  described  in  Chapter  V,  a  low  resolution 
simulation  was  used  to  forecast  the  component  failures  associated  with  the  Weibull 
renewal  process.  The  advantage  of  such  simulations  is  the  ability  to  include  other 
factors  which  may  afiect  a  component's  failure  rate.  Examples  of  these  factors  would 
include  number  and  type  of  landings,  mission  type,  maintenance  programs,  weather 
conditions,  etc.  It  is  safe  to  say.  that  these  more  coniplicated,  high  resolution  models 
could  include  more  information  about  the  factors  which  cause  components  to  fail,  and 
should  ultimatelv  result  in  better  forecasts  of  demand. 
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APPENDIX  A 
TRANSACTION  CODES 

CODE    DESCRIPTION 

00  Used  to  report  an  inventor>'  gain. 

02  Used  to  report  a  change  in  the  readiness  reportable 
status  of  an  equipment. 

03  Used  to  report  an  equipment  loss. 

11  a.  On-equipment  work  not  involving  a  removal  or 
replacement  of  a  defective  component  or  item. 

b.  On  supprotine  ensine  documents  not  havins  a  removal 
of  a  defective  or^suspected  defective  component  or  item 
when  the  ensine  is  not  specifically  identified  to  a 

•  specific  aircraft. 

c.  Used  at  the  O  or  I  level  maintenance  activities  when 
closing  out  a  maintenance  action. 

12  On-equipment  work,  including  eneines.  involving  non- 
repairable  components  items  aocifmented  as  failed  parts. 

14  Removal  of  a  nondefectiye  component  item,  excluoins 
cannibalization.  from-an  eneine  to  be  processed  at  the 
O-level  maintenance  activity. 

15  Installation  of  a  nondefective.  component,  item. 
eJccludine  cannibalization.  on  an  engine  to  be  processed 
at  an  O-Tevel  maintenance  activity.  " 

16  Removal  of  a  nondefective  component/item,  excluding 
engine  components  items  and  a  cannibalization  to  be 
processed  at  an  O-level  maintenance  activity. 

17  Installation  of  a  nondefective  component/item 
excluding  engine  component  items  and  cannibalization. 

18  a.  The  removal  and  replacement  of  nondefective 
engines  and  component,  items  to  accomplish  a 
cannibalization  action. 

b.  The  removal  and  replacement  of  those  consumeable 
components  items  subject  to  a  scheduled  removal  interval 
or  items  of  supplv  significance. 

c.  The  removal  and  replacement  of  a  nondefective  engine 
component  for  cannibalization  at  the  1-level  only. 

19  a.  The  removal  and  replacement  of  nondefective 
engine  component  to  accomplish  a  cannibalization  at  the 
O-Tevel  onlv. 

b.  The  removal  and  replacement  of  those  consumeable 
engine  components  sumect  to  a  scheduled  removal  interval 
or  Items  of  supply  significance. 

21  Will  be  used  when  a  repairable  component  is  removed. 

excluding  engine  components,  for  processing  at  an  I  or  D 
level  niafntenance  activitv  (This  code  is  used  when 
only  the  removal  must  be  documented  and  a  replacement 
is  not  required). 

23  Removal  and  replacement  of  a  defective  or  suspected 

defective  component  from  an  end  item,  excluding  engines 

at  the  O-lcvel. 

Additionally,  this  code  is  used  for  the  removal  and 
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replacement  of  a  complete  engine  assembly  for  a  defect, 
suspected  defect,  or  scheduled  maintenance  requirement. 

24  Will  be  used  when  a  repairable  engine  component 
is  removed  for  processing  at  an  I  or  D  level  activity 
(This  code  is  used  when  onlv  the  removal  must  be 
documented  and  a  replacement  is  not  required). 

25  Removal  and  replacement  of  a  defective  or  suspected 
defective  repairaole  component  from  an  engine. 

30  Used  to  document  components  processed  through  the 
I-level  maintenance  activity  for  check,  test,  and" 
service. 

31  Work  performed  on  a  removed  repairable  component  with 
no  failed  parts  or  awaitine  parts  documented  in  the 
Failed/Required  Material  mocks  of  the  VIDSMAF. 

32  Work  performed  on  a  removed  repairable  component  with 
failed  parts,  awaitine  parts,  or  cannibalization 

actions  documented  m  the  Failed/Required  Material 
blocks. 

39  Close  out  for  man  hours  or  awaitine  parts  at  the 

I-level  maintenance  activity. 

41         .  Technical  directive  compUance  with  no  part  number 
change. 

47  Technical  directive  compliance  with  a  part  number 

change. 

72  Will  be  useed  to  report  SCIR  data  bv  the  reportine 

custodian  when  transient  maintenance  is  peribrmea  bv 
other  than  the  reporting  custodian. 
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APPENDIX  B 
ACTION  TAKEN  CODES 

CODE    DESCRIPTION 

1  BCM  -  maintenance  activity  not  authorized  to  conduct 
repair 

2  BCM  -  Lack  of  equipment,  tools,  or  facilities 

3  BCM  -  Lack  of  technical  skills  necessary  to  complete 
repair 

4  BCM  -  Lack  of  parts 

5  BCVI  -  Fails  check  and  test,  and  maintenance  is  allowed 
to  conduct  check  and  test  only 

6  BCM  -  Lack  of  technical  data 

7  BCM  -  Beyond  the  authorized  repair  depth  ; 

8  BCM  -  Administrative 

9  BCM -Condemned,  repair  not  feasible 

A  Items  of  repairable  Vlaterial  or  Weapon  Support 

Svstem  Discrepencv  Checked  No  Repair  Required.  This 
code  IS  used  for  ail'discrepencies  which  are  checked 
and  found  that  either  the  reported  deficiencv  cannot 
be  duplicated,  or  is  operating  within  allowable 
tolerances. 

B  Repair  or  replacement  of  attaching  units,  seals.- 

gasKets.  etc..  that  are  not  integral  parts  of  work 
unit  coded  items. 

C  This  code  is  used  when  a  repairable  item  of  material 

identified  by  a  work  unit  code  is  repaired. 

D  This  code  is  used  to  closeout  a  VIDSMAF  when 

component  repair  is  to  be  peribrmed  at  another  facility. 

F  Failure  of  components'items  undergoing  test. 

J  This  code  is  used  when  an  item  is  calibrated  and  found 

serviceable  without  need  for  adjustment. 

K  This  code  is  used  when  an  item  must  be  adjusted  to 

meet  calibration  standards. 

L  This  code  is  used  when  a  maintenance  action  must 

be  stopped  or  delaved  while  awaiting  parts  which  are 
not  available  locallv.  and  a  component  goes  into 
an  awaiting  parts  status. 

N  This  code  is  used  bv  an  organizational  activitv 

when  it  becomes  ne"cessar\^to  closeout  a  maintenance 
action  during  or  at  the  end  of  a  reporting  period 
for  any  reason. 

P  This  code  is  entered  when  an  item  of  material  is 

removed  and  onlv  the  removal  is  to  be  accounted  for. 
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Q  This  code  is  entered  when  an  item  is  installed  and 

only  the  installation  action  is  to  be  accounted  for. 

R  This  code  is  entered  when  an  item  of  material  is 

removed  due  to  suspected  malfunction  and  the  same  or  a 
like  item  is  reinstalled. 

S  This  code  is  entered  when  an  item  of  material  is 

removed  to  facilitate  other  maintenance  and  the  same 
item  is  reinstalled. 

T  This  code  is  used  when  an  item  of  material  is 

removed  and  replaced  for  a  cannibalization  action. 

Y  This  code  is  used  when  the  time  expended  in  locating 

a  discrepency  is  great  enough  to  warrant  separating 
the  troublesnoot  time  from  the  repair  time. 

Z  This  code  is  used  when  actually  treating  corroded 

items,  and  includes  cleaning,  treatment^  priming  and 
paintmg. 
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APPENDIX  C 
WHEN  DISCOVERED  CODES 

CODE      DESCRIPTION 

A  Before  flight  -  mission  aborted  -  failure  discovered  by  aircrew 

B  Before  flight  -  no  abort  -  failure  discovered  by  aircrew 

C  Inflight  -  mission  aborted 

D.  Inflight  -  no  abort 

E  After  flight  between  flight  -  failure  discovered  by  aircrew 

F  Pilot/NFO  weekly  inspection 

G  Acceptance/transfer  inspection 

H  Between  flights  -  failure  discovered  by  ground  crew 

J  Daily  inspection 

K  Preflight.  daily; preflight,  postflight,  or  turnaround  inspection 

L  Special  inspection 

M  Calender  odd/major/phase  inspection 

N  Calender  even  inspection 

O  Administrative 

P  Functional  checkflight 

Q  Conditional  inspection 

R  Quality  assurance  inspection 

S  Oil  analysis  inspection 

U  Modification,  SDLM,  Overhaul/ Airline  maintenance 

V  Related  maintenance  action 

W  In-shop  repair; disassembly  for  maintenance 

X  Test  bench/engine  test  stand  operation 

Y  Receipt  or  withdrawal  from  supply 
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APPENDIX  D 
MAXIMUM  LIKELIHOOD  DERIVATIONS 

1.        MAXIMUM  LIKELIHOOD  ESTIMATORS 

a.  Properties 

The  method  used  in  this  thesis  to  estimate  the  model  parameters  is  maximum 
likelihood.  As  such,  it  is  necessary  to  review  two  important  properties  of  this  estimator. 
First,  if  a  random  sample  of  size  n  is  taken  of  a  random  variable  X  whose  probability 
mass  function,  p{x),  depends  on  an  unknown  parameter,  0,  then  the  probability  law  for 
the  maximum  likelihood  estimator  is  asymptotically  normal  (as  n  gets  large)  with: 

MEAN  =  e   and  VARIANCE  =  (nK^)-^  (eqn  D.l) 

where,  K'^=  E[{(5/^G)lnp(x))^]  for  a  discrete  random  variable  X,  and  n  equals  the  total 

number  of  observations.    Thus,  the  maximum  likelihood-  estimator  is  asymptotically 

unbiased:   that  is,   the   expected   value   of  the   estimator  is  equal   to   the   unknown 

parameter  0.  Furthermore,  the  variance  of  the  estimator  is  the  Cramer-Rao  lower 

bound,    and    represents    the    smallest    possible    variance    of  an   unbiased    estimator. 

[Ref  10:  pp.  372,379] 

The  second  property  describes  the  probability  law  for  maximum  likelihood 

estimators  when  the  X.'s  for  i=  1 n  are  not  identicallv  distributed.  In  the  models 

1  ' 

discussed  in  this  thesis,  the  probability  mass  function  for  each  X.  was  dependent  on 
one  or  more  other  random  variables.  Nevertheless,  the  distribution  of  the  maximum 
likelihood  estimator(s)  is  still  asymptotically  normal  with: 

MEAN  =  0   and  VARIANCE  =  (Y Kj^)'!  (eqn  D.2) 

i 

where,  Kj'-=  E[((^/^0)lnp(x))"]  is  evaluated  using  the  specific  parameters  defining  the 
i  observation.  The  difTerence  in  the  variance  in  equation  D.l  and  the  variance  in 
equation  D.2  is  caused  by  probability  mass  functions  for  the  X-'s  which  are  not 
identically  distributed.  For  a  complete  derivation  of  this  result,  the  reader  is  referred  to 
W.  J.  Heintzelman's  1975  paper  on  Determining  the  Failure  Rate  When  Failure  Times 
Are  Not  Known  Exactly  [Ref  llj. 
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b.  Confidence  Intervals 

Confidence  intervals  for  the  unknown  parameters  can  then  be  calculated  using 
the  maximum  likelihood's  probability  law  described  in  the  previous  section.  These 
limits  are  defined  as: 

CONFIDENCE  LIMITS  =  6'±(zi.a/2  ^  ^6  ) 

where,  8'  is  the  maximum  likelihood  estimator  for  0,  Cq'  is  the  standard  deviation  of 
the  maximum  likelihood  estimator,  and  Zj_q/2  is  the  l-a/2  quantile  of  the  standard 
normal  distribution. 

c.  Formulation  of  Probability  Mass  Functions 

The  generalized  format  used  to  construct  each  model's  probability  mass  function  for  X- 
is: 

p(.x.)  =  [P{ failure  occurred)]    '^i  x  [P( right  hand  censored)]^! 

The  random  variable,  X.,  is  an  indicator  variable  equal  to  zero  if  the  data  entr\' 
contains  a  failure,  and  equal  to  one  if  the  observation  is  right  b.and  censored.^  The 
following  variables  are  used  to  derive  the  maximum  likelihood  estimate(s)  in  each  of 
the  models: 


•  n 

•  m 

•  n-m 


Total  number  of  records  in  the  final  data  set 

Total  number  of  records  containing  component  failures 

Tot^l  number  of  records  which  were  right  hand  censored 


2.        EXPONENTIAL  MODEL 

a.  Properties 

The  probability  density  function  f(t),  survival  function  S(t),  and  hazard  rate 
function  r(t)  for  the  Exponential  distribution  are  defined  below. 


S(t)=Pr(T>t)=e*^^ 


r(t)  =  ililx 
S(t} 


^Risht  hand  censored  observations  identify  those  components  that  did  not  fail 
to  the  termination  of  the  deployment. 


prior 
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b.  Formulation  of  the  Maximum  Likelihood  Equation 

This  section  outlines  the  formulation  of  the  maximum  likelihood  equation 
used  to  find  X',  the  maximum  likelihood  estimate  for  the  failure  rate  parameter  X.  The 
final  data  set  contains  information  supporting  two  record  types.  The  first  record  type 
identifies  observations  of  components  failing  during  some  specified  interval  of  time.  T- 
to  T-  +  t-.  The  probability  equation  representative  of  this  type  of  record  is  given  in 
equation  D.3. 

P(Tj<t<Ti  +  tj)=exp(-XTi)x(l.exp(-Xti))  (eqn  D.3) 

The  second  record  type  pertains  to  an  observation  of  a  component  surviving 
past  some  time  T-.  These  censored  observations  contain  no  information  regarding  the 
actual  time  of  component  failure.  Equation  D.4  defines  the  probability  relationship 
associated  with  this  type  of  observation. 

■S(T.)=exp(-XTj)  .  (eqn  D.4) 

The  probability  mass  function,  p{Xj),  for  the  Exponential  model  is  obtained  by 
combining  equations  D.3  .and  D.4.  The  random  variable  X.  equals  zero  when  a  failure 
is  observed  and  one  when  the  data  is  right  hand  censored  (i.e.,  component  survived 
until  the  end  of  the  deployment). 

p(Xi)  =  [exp(-XT.)r^i  X  [exp(-XTi)  x  (l.exp(-Xtj)/l-V  (eqn  D.5) 

where. 


Xj=[0,l},  X>0,  andTj.  tj>0 

T-  equals  the  aggregate  flight  time  \vithout  component  failure  for  the  ith  record 

tj  equals  the  fiight  duration  of  the  sortie  during  which  the  failure  was  observed 


The  likelihood  function  for  p(Xj),  L(p(x-))  can  be  expressed  as  the  product  of 
the  n  probability  mass  functions.  The  maximum  likelihood  estimate  for  X  is  the  value 
of  X'  which  maximizes  L(p(x-)).  Since  the  logarithmic  function  is  a  monotonically 
increasing  function,  the  X'  which  maximizes  L(p(x-))  is  identical  to  the  X'  which 
maximizes  the  log-likelihood  function  K(p(x-)).  The  optimal  value  for  X'  is  found  by 


taking  the  derivative  of  the  log-likelihood  function  with  respect  to  X,  setting  it  equal  to 
zero,  and  solving  for  X. 

Likelihood  function: 

L(p(Xj))  =  n  [exp(->.Tj)]'^i  X  [exp(-XTi)  x  (l.exp{-Xtj)](l-^i)      (eqn  D.6) 

i=  1,2,3, .••,n 

Log-likelihood  function: 

K(p(Xi))  =  X  -  XJ.x..  XT.{\-x.)  +  (l-Xj)log[l-exp(.Xtj)]  (eqn  D.7) 

i 

The  partial  derivative  of  the  log-likelihood  function  with  respect  to  X  is: 

dK    ^  [(l-x.H.exp(-Xt.)] 

_  =  y  .T.x. -T.(l-x.)  +  ^^^ '-^-^ — '—  (eqnD.8) 

^X     f       •  •       •'      ''  [l-exp(-Xtj)]  ^ 

Equation  D.9  is  the  result  of  setting  Equation  D.8  equal  to  zero  and 
rearranging  the  summations  to  eliminate  the  x-(s).  The  left  hand  side  of  Equation  D.9 
is  summed  over  the  set  (i)  of  all  records,  censored  and  uncensored.  But,  the  right  hand 
side  is  only  summed  for  the  uncensored  observations  (j). 

^         ^  t.exp(-Xt.)  T;,t:>0.  X>0,  m<n 

Vt.  =  V_lJ_! L_  1']  (eqn  D.9) 

Y     '     Y  l-exp(-Xt)  1=1.2 n    j=l,2,...,m 

c.  Maximum  Likelihood  Variance 

As  stated  in  Section  D.l,  the  maximum  likelihood  estimator  for  X  is 
approximately  normally  distributed  with  mean  equal  to  X  and  variance  equal  to  1  '^K." 
for  i=  l,...,n.  In  this  model: 

(t.exp(-Xt.)  -  T.)^  -y 

Kf  = ■  ,        ;, ' X  [exp(-XTi)  X  exp(  -Xtp]  +  T-exp(-XTi) 

'  l-exp(-Xt')  ^  1         i  1 

where,  Kj^  is  evaluated  using  the  maximum  likelihood  estimate  of  X  (i.e.,  X  =  X').  The 
standard  deviation  of  the  maximum  Ukelihood  estimator  and  associated  confidence 
limits  displayed  in  Table  7  were  computed  using  the  equation  listed  above. 
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3.       WEIBULL  MODEL 

a.  Properties 

The  probability  density  function  f^t),  survival  function  S(t),  and  hazard  rate 
function  r(t)  for  the  Weibull  distribution  are  displayed  below. 

n;t)=>.*^at^-^exp[-(Xt)"]         X,  a>0    t>0 

S(t)  =  exp[-{>.t)«] 

r(t)  =  X«at«-^ 

b.  Formulation  of  the  Maximum  Likelihood  Equation 

This  section  describes  the  formulation  of  the  equation  used  to  find  the 
maximum  likelihood  estimates  of  the  Weibull  scale  (X)  and  shape  (a)  parameters.  The 
two  different  record  types  supporting  this  equation  are  identical  to  those  identified  in 
the  preceding  section.  The  key  difference  between  the  two  models  is  that  the  Weibull 
distribution  possesses  the  "memor\'less"  property  for  a  equal  to  one  only.  As  such,  the 
conditional  survival  probability  is  usually  not  equal  to  the  unconditional  survival 
probability.  The  conditional  survival  probability  for  the  Weibull  distribution  is  defined 
asi^o         '  . 

^^^_exp{-[X(T.  +  t.,)]«] 
't/V-      exp[-(XT.)«] 

The  conditional  probability  of  failure  during  the  interval  T.  to  T.  + 1.  can  be 

expressed  as  l-Sj(tj),  the  complement  of  the  conditional  survival  probability.    The 

i 
probability  mass  function  for  the  Weibull  model  is  then  formed  by  combining  the 

conditional  probability  of  a  component  failure  within  the  interval  Tj  to  Tj-i-t-  (x-  =  0) 

with  the  probability  that  a  component  survives  time  T-  (x- =  1). 

p{Xi)=[l-S-j-_(ti)]l-\  X  exp[-{XTj)Yi  (eqn  D.ll) 

where, 


• 


x.=  {0,l},  X>0,  andT..  t.>0 

Tj  equals  the  aggregate  fiight  time  without  component  failure  for  the  ith  record 


.    The   subscripted  T:   is   used   to    indicate   the   conditional   survival   probabilitv 
associated  with  a  component  that  has  survived  T:  hours  of  operation. 
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•      t-  equals  the  flight  duration  of  the  sortie  during  which  the  failure  was  observed 
The  log-likelihood  function,  K(p(Xj)),  is: 

K(p(Xj))  =  y  (l-x.)(log[l-S-p{tj)]}  +Xi[-(XT/]  (eqn  D.12) 

i  * 

i=  l,2,3,...,n 

The  partial  derivatives  of  the  log-likelihood  function,  K(p(x.)),  with  respect  to  the 
parameters  X  and  a  are  displayed  below.  They  have  been  simplified  by  removing  the 
x-(s)  and  summing  over  the  appropriate  records.  The  set  i=  1,2, 3,. ..,n  represents  all  data 
records  (censored  and  uncensored),  and  the  set  j  =  m+  l,m+2,...,n  includes  only  those 
records  which  are  censored. 

aK  ^       aX«-'|(T,  +  t,)«-Ti«|ST.(t,)      ^^^„.,^„ 

dx   7  i-Sjitj)  Y         ' 

i  ■•  =  . 

eK_^     i[X-(T,  +  ti)|«loglX(Ti  +  ti)l-(>.T,)"log(?.T,)S-r.(t,)        V(),T)«lo«).T) 

.sa    Y  i-s-r.it;)     ■,  f   y    -^   y 

c.  Maximum  Likelihood  Variance 

The  following  equations  were  used  to  compute  the  variance  for  the  maximum 
likelihood  estimators,  X,'  and  a'.  These  variances  were  then  used  to  construct  the 
confidence  limits  displayed  in  Table  8. 

K-^  with  respect  to  X': 

Ki2(  V)  =  («X«- '  ,2  X  {    '"Ti^V"-T."'ST(V!-        ^  T>exp[.(XTi)«l  j 

i 

K-    ■with  respect  to  a': 

.         ([[MTi  +  ti)]«log[X(Ti  +  ti)HXTi)«log(XTi)}Sy(tj))2 
K:"(a  )  = — — *-i— ^ p 

i 

where,  p  equals  [(XT|)"log(XT-)]^exp[-{XT^)"j. 
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4.        GEOMETRIC  MODEL 

a.  Properties 

The  probability  mass  function,  Pr{F=s),  and  survival  function,  Pr(F>s),  for  the 
Geometric  distribution  are  given  below. 

Pr(F=s)  =  pq^"^       0<p<l,  q=l-pand  s=  1,2,3,.. . 
Pr(F>s)=q^'^ 

where,  F  is  the  random  variable  representing  the  observation  of  the  first  failure,  and  q 
is  the  probability  of  the  component  surviving  the  sortie.  Equations  D.13  and  D.14 
illustrate  the  Geometric  distribution's  "memor>'less"  and  constant  failure  rate  properties 
respectively. 

Pr{F>a4-b/F>a)  =  (q^^^Vq''^  =  q^*  =  Pr,(F>b)  (eqn  D.13) 


Pr(F=s);Pr(F>s)  =  (pq^-^)/q^"^  =  p  (eqn  D.14) 


b.  Formulation  of  the  Maximum  Likelihood' Equation 

This  section  formulates  the  equation  used  to  find  the  maximum  likelihood 
estimate  of  the  constant  failure  probability  (p)  associated  with  any  given  sortie.  For 
the  Geometric  model,  the  two  event  types  which  can  occur  include  components  failing 
during  the  s  sortie,  or  components  failing  during  some  unknown  sortie  after  the  s^ 
sortie.  The  probability  expression  associated  with  the  first  type  of  event  is  given  in 
equation  D.15,  and  equation  D.16  defines  the  relationship  for  those  observations  where 
the  component  survived  at  least  s  sorties. 


Pr(F=s)  =  pq^-1  (eqn  D.15) 


—  ^s 


Pr(F>s)  =  q^  {eqn  D.16) 

The  probability  mass  function  for  this  model  is  obtained  by  combining 
equations  D.15  and  D.16.  Again,  X  is  allowed  to  take  the  values  zero  and  one  only, 
with  zero  indicating  a  component  failure  on  the  s  sortie  and  a  one  indicating  a 
component  failure  after  the  s      sortie. 
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p(Xi)=[p(l-p)VA]A-^  X  [{\.p)h]\  (eqn  D.17) 

where, 

•  0<p<l,q=l-p,  s-=  1,2,3,... 

•  s-  equals  the  aggregate  number  of  sorties  flown  inclusive  of  the  sortie  causine 
the  failure  when"  xp=  0  and  equals  the  aggregate  number  of  sorties  without  a 
component  failure  when  x-=  1 

The  maximum  likelihood  estimate  p'  of  the  unknown  sortie  failure  rate  will  be 
found  by  taking  the  derivative  of  the  log-likelihood  function  [K(p(x-))]  with  respect  to 
p,  setting  the  derivative  equal  to  zero  and  solving  for  p. 

K{p(Xi))  =  V  (1-Xi)ln{p)  +  {Si-l)(I-xpin(l-p)  +  s-Xjln(l-p)       (eqn  D.18) 
i 

i=  1,2.3 n 

^K_     (l-x^)     [(s--l)(l-x-)]        (S-X-) 


.      dp     Y    P  <1-P>  (1-P^ 

The  partial  derivative  listed  above  can  be  simplified  by  removing  the  indicator 
variables  and  adjusting  the  suminations  to  produce  the  closed  form  expression  for  p' 
given  in  equation  D.19. 

vi  +  _l.  =  v_!L 
Y  p      i-p     T  ^"P 


m 


In 


i=l,2,...,n  (eqn  DT9) 


where,  j=  1,2,3,. ..,m  is  the  set  of  uncensored  records  (i.e.,  component  failures),  and 

i=  1,2,3 n  is  the  set  of  all  records,  censored  and  uncensored.    Thus,  the  maximum 

likelihood  estimator  p'  equals  the  number  of  flights  experiencing  a  component  failure 
divided  by  the  total  number  of  sorties  flown, 
c.  Maximum  Likelihood  Variance 

The    maximum    likelihood    estimator's    variance    was    computed    using    the 


formulas   described   in   Section   DT.    The   equation   defining   the   Geometric   model's 
i 
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expression  for  K/-  is  given  below: 


Ki^  =  [(l-psi)2(l.p)V3]/p  +  Sj2(l.p)V2  (eqn  D.20) 

where,  K-  is  evaluated  at  p  =  p'.  These  variances  were  used  to  compute  the  confidence 
limits  listed  in  Table  9. 

5.        MIXED  MODEL 

a.  Properties 

Combining  a  constant  sortie  failure  probability,  p^,  with  a  flight  hour  probability 
distribution  results  in  an  expression  for  the  probability  of  failure,  p-,  for  the  i^^  sortie. 
This  is  expressed  as  the  complement  of  a  component  surviving  both  the  stresses  of  the 
sortie  and  the  flight  hours  flown: 


p-=  l-exp(->.t-)(l-pQ)       0<Po<l.X>0,t^>0  (eqnD.21) 

where,  the  flight  hour  failure  distribution  is  Exponential. 

The  model  is  then  formulated  in  the  same  manner  as  the  Geometric  model 


describe4  in  the  preceding  section.  The  probability  that  the  first  failure  occurs  on  the 
s      sortie  is  e; 
for  i=  l,...,s-l. 


s      sortie  is  expressed  in  equation  D.22.    Note  that  the  p^^^.  are  not  necessarily  equal 


Pr(F=  s)  =  p_  n  (I-pi)        i=  l,2,...,s-I  (eqn  D.22) 

i 


Equation  D.23  results  from  using  the  expression  for  p-  given  in  equation  D.21 
to  rewrite  equation  D.22. 

Pr(F=s)  =  [l-exp(-Xt)(l-pQ)][(l-pQ)S-lexp(-XT)]       T-,t->0     (eqn  D.23) 

b.  Formulation  of  the  Maximum  Likelihood  Equation 

This  section  describes  the  formulation  of  the  equation  used  to  find  the 
maximum  likelihood  estimates  for  the  failure  rate  parameter  (X)  and  the  constant  sortie 
probability  of  failure  (p^).  The  probability  relationship  supporting  the  observation  of 
the  first  failure  on  the  s*-"  sortie  is  given  in  equation  D.23,  and  the  expression  for  the 
observation  of  a  failure  at  some  unknown  time  after  s      sortie  is  listed  below. 
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Pr{F>  s)  =  n  (1-p:)  =  exp(-X.T){l.p_)S  (eqn  D.24) 

1  '■  ^ 


The  probability  mass  function  p(Xj)  for  this  model  results  from  combining 
equations  D.23  and  D.24.  Notice  that  when  X  equals  zero  the  probability  mass 
function  for  the  Mixed  model  is  identical  to  the  probability  mass  function  expressed  for 
the  Geometric  model  in  equation  D.17,  and  when  p^  equals  zero  the  probability  mass 
function  is  equivalent  to  the. Exponential  probability  mass  function  given  in  equation 
D.5.  The  above  statement  implies  that  the  Mixed  model  has  the  capability  of 
modelling  the  Geometric  or  Exponential  cases  as  the  situation  would  require. 

p{Xj)=  [[l-exp(-Xti)(l-pQ)][(l-pQ)Srlexp(.XT^)]}l-^  x  [(l.p^)^-,  exp(-XT^)]^ 

where. 

•  0<Pq<  1,  ;.>0,  s—  1,2,3 and  J.,  tj>0 

•  T-  equals  the  aggregate  flight  time  without  component  failure  for  the  ith  record 

•  tj  equals  the  flight  duration  of  the  sortie  during  which  the  failure  was  observed 

•  S;  equals  the  aggregate  number  of  sorties  flowTi  inclusive  of  the  sortie  causing 
the  failure  when  x:'"=  0  and  equals  the  aggregate  number  of  sorties  without  a 
component  failure  when  x-=  1. 

The  log-likelihood  function  K(p(x-)),  and  the  partial  derivatives  of  the  log- 
likelihood  functions  with  respect  to  the  parameters  X  and  p^  are  given  below.  The 
maximum  likelihood  estimates  for  "k  and  p^  were  obtained  utilizing  the  Quasi-Newton 
method  described  in  Appendix  G.    The  log-likelihood  function,  K(p(x-)),  is: 


K(p(x-))  =  y  ( l-xp{ln[  l-exp(-Xt-)(  I-Pq)]  +  [(s--l)ln(l-Po)]-XTi}  +  x- {[Sjln(l-Po)]-XT^} 
i 
i=l,2,3,...,n 


The  partial  derivatives  of  the  log-likelihood  function,  K(p(x.)),  have  been  simplified  by 

removing  the  x-(s)  and  summing  over  the  appropriate  records.  The  set  i=  1,2,3 n 

represents  all  data  records  (censored  and  uncensored),  and  the  set  j  =  m+ l,m-f- 2,...,n 
includes  only  those  records  which  are  censored. 
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dK          (l-pQ)t-exp(-Xt-)      _     ^  _  /         .  o.^ 

TT'^zl  -, ,  ;    ,,,     \  -  Ti  -  X  Tj  (eqn  A.25) 

^K        ^  eXp(-Xt:)  S;-l  ^  S: 

—  =  y — ^- E ^—  (eqn  A.26) 

^Po      i   l-exp(->.tp(l-pQ)       1-Pq        j      1-Pq 


c.  Maximum  Likelihood  Variance 

The  equation  for  the  variance  of  the  probability  distribution  for  p^'  is 
identical  to  equation  D.20  when  X  equals  zero.  The  equation  for  the  K-  for  "k'  is  given 
below. 

Ki2(X')  =  [{{l-p^)t^]/PQ.Ti}2xp^(l.p^)Si-l  +  Ti2(l-p^)Si 


7 

where,   K-"  is  evaluated  at  'k  =  \'  =  0.    These  variances  were  used  to  compute  the 
confidence  limits  for  X  listed  in  Table  9. 
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APPENDIX  E 
SINGLE  PARAMETER  OPTIMIZER 

*  COMPUTES  MAXIMUM  LIKELIHOOD  VALUE  FOR  A  SINGLE  UNKNOWN 

*  PARAMETER  USING  GOLDEN  RATIO  LINE  SEARCH. 

IMPLICIT  REAL*8  (A-H,0-Z) 
COMMON /RL/FTNF , FTF , PARAM , PSD 
COMMON/ INT/CENSOR , NF , COUNT 

DIMENSION  CENSOR(200) ,NF{200) , FTNF(200) ,FTF(200) ,PARAM(2) 
INTEGER  CENSOR, COUNT, TFLAG 
•REALMS  LB,LC 

DATA  LB/.OOOlDO/ITER/1/K/l/STEP/.OOOOlDO/ALPHA/.OOOlDO/ 

*  INITIALIZE 

C0UNT=1 
PARAM(1)=.00001D0 

TFLAG=0 

*  COMPUTATION  OF  SEARCH  INTERVAL  TAU.  TAU  IS  APPROXIMATELY 

*  .618  AND  TAUl  IS  APPROXIMTELY  .382 

TAU=(5*'^.5-l)/2 

TAU1=1-TAU  ■  ■  .   . 

*  INPUT  DATA 

5  READ  (3,50,END=99)  CENSOR ( COUNT ) ,NF( COUNT ), FTNF( COUNT ), FTF (COUNT) 

C0UNT=C0UNT+1 

GO  TO  5'    "  .   • 

99     C0UNT=C0UNT-1 

*  BOUNDING  PHASE  -  SWANN'S  METHOD 

CALL  ML(LB,FI) 

6  ALPHA=ALPHA+(2'^*K)*STEP 

CALL  ML ( ALPHA, F2) 

*  IF  NEW  VALUE  IS  LESS  THAN  PREVIOUS  VALUE  TERMINATE  BOUNDING 

IF  (F2.LT.F1)  GO  TO  10 
F1=F2 
K=K+1 
GO  TO  6 
10     UB=ALPHA 

^DEFINE  INITIAL  POINTS  &  ASSOCIATED  FUNCTIONAL  VALUES 

B1=TAU1*(UB-LB)+LB 
B2=TAU*(UB-LB)+LB 
CALL  ML(B1,F1) 
CALL  ML(B2,F2) 

*  COMPARES  THE  FUNCTIONAL  VALUES  OF  THE  INTERNAL  PTS .  Bl  &  B2 ,  AND 

*  SELECTS  THE  VALUE  WITH  THE  LARGEST  FUNCTIONAL  VALUE  AS  THE  UPDATED 

*  PARAMETER  ESTIMATE 

150    IF  (F1.GT.F2)  THEN 
PARAM(2)=B1 
UB=B2 
B2=B1 
F2=F1 

B1=TAU1*(UB-LB)+LB 
CALL  HL(B1,F1) 
ITER=ITER+1 
GO  TO  200 
END  IF 
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IF  (F2.GE.F1)  THEN 

PARAM(2)=B2 

LB=B1 

B1=B2 

F1=F2 

B2=TAU*(UB-LB)+LB 

CALL  ML(B2,F2) 

ITER=ITER+1 

GO  TO  200 
END  IF 

200    CALL  TERM(F1,F2,TFLAG) 
PARAM(1)=PARAM(2) 
IF  (TFLAG.EQ.l)  GO  TO  275 
GO  TO  150 

*   WRITE  OUTPUT 

275   WRITE  (6,300) 

WRITE  (6,310)  PARAM(l) 

CALL  SD 

WRITE  (6,305)  PSD 

UC=PARAM ( 1 ) +1 • 96*PSD 

LC=PARAM(1)-1.96*PSD 

WRITE  (6,320)  LC,UC 


* 

FORMAT 

50 

FORMAT 

300 

FORMAT 

305 

FORMAT 

310 

FORMAT 

320 

FORMAT 

9999 

STOP 

END 

'I1,1X,I4,5X,F6.1,1X,F3.1) 

■ 'OWUC: 7'OCONVERGENCE  CRITERION  ESTABLISHED') 

'OSTANDARD  DEVIATION  IS        :  ',F9.8) 

'OMLE  FOR  UNKNOWN  PARAMETER  IS:  ',F9.8) 

'0.95' ASSYMPTOTIC  CONFIDENCE  LIMITS  ARE  : ' , F5 .4 , ' , ' , F5 .4) 


SUBROUTINE  TERM(F1 , F2 ,TFLAG) 

DETERMINES  IF  TERMINATION  CONDITION  HAS  BEEN  SATISFIED 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/FTNF , FTF , PARAM , PSD 

COMMON/  INT/ CENSO^^  ,  NF ,  COUNT 

DIMENSION  CENSOR(200) ,NF(200) , FTNF(200) , FTF(200) ,PARAM(2) 

INTEGER  TFLAG 

TEST=ABS(F1-F2) 

IF  (TEST. LE. .00003)  TFLAG=1 

RETURN 
END 

SUBROUTINE  ML(PT,FVAL) 

COMPUTES  LOG-LIKELIHOOD  FUNCTIONAL  VALUE 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/FTNF , FTF , PARAM , PSD 

COMMON/ INT/ CENSOR , NF , COUNT 

DIMENSION  CENSOR(200) ,NF(200) , FTNF(200) ,FTF(200) ,PARAM(2) 

INTEGER  CENSOR, COUNT, TFLAG 

REALMS  Kl 

FVAL=0 

INPUT  FLIGHT  DATA 
DO  5  1=1, COUNT 

UNCENSORED  RECORDS 
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IF    (CENSOR(I).EQ.O)    THEN 

K1=EXP(-PT*FTF(I)) 

FVAL=FVAL+LOG ( 1 -Kl ) -PT*FTNF ( I ) 

GO   TO    5 
END    IF 

*  CENSORED  RECORDS 

FVAL=FVAL-PT*FTNF ( I ) 

5  CONTINUE 

RETURN 
END 

SUBROUTINE  SD     '         " 

*  COMPUTES  ASSYMPTOTIC  STANDARD  DEVIATION  FOR  THE  MAXIMUM  LIKELIHOOD 

*  ESTIMATE 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/ FTNF , FTF , PARAM , PSD 

COMMON/ INT/ CENSOR , NF , COUNT 

DIMENSION  CENSOR(200) ,NF(200) , FTNF(200) ,FTF(200) ,PARAM(2) 

INTEGER  CENSOR, COUNT 
REALMS  KISQ,K1,K2,K3 

SUMKI=0 . 

*  INPUT  FLIGHT  DATA 

DO  5  1=1, COUNT 

*  UNCENSORED  RECORDS 

IF  (CENSOR(I).EQ.O)  THEN 

K1=EXP ( -PARAM ( 1 ) *FTF ( I ) ) 

K2=EXP ( -PARAH ( 1 ) *FTNF ( I) ) 

K3=((FTF(I)^K1/(1-K1))-FTNF(I))*((FTF(I)'^K1/(1-K1))-FTNF(I)) 

KI S0=K3^  ( K2*  ( 1  -Kl ) )  +  ( FTNF  ( I )  =^FTNF  ( I )  *K2 ) 

GO  TO  6 
END  IF 

*  CENSORED  RECORDS 

KISQ=FTNF  ( I )  '*^FTNF  ( I )  *K2 

6  SUMKI=SUMKI+KISQ 
5      CONTINUE 

PSD=SQRT(1/SUMKI) 

RETURN 
END 

SUBROUTINE  SDl 

*  COMPUTES  ASSYMPTOTIC  STANDARD  DEVIATION  FOR  THE  MAXIMUM  LIKELIHOOD 

*  ESTIMATE  FOR  THE  SORTIE  DEPENDENT  MODELS 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/ RL/ FTNF , FTF , PARAM , PSD 

COMMON/ INT/ CENSOR , NF , COUNT 

DIMENSION  CENSOR(200) ,NF(200) ,FTNF(200) ,FTF(200) ,PARAM(2) 

INTEGER  CENSOR, COUNT 

REALMS  L,LL1,LL2 

*  INITIALIZE  PARAMETER  ESTIMATES 

P=.0232D0 
L=.OOOODO 
SUMKI1=0.0D0 
SUHKI2=0.0D0 
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*  INPUT  FLIGHT  DATA 

DO  5  1=1, COUNT 
C1=(1-P*NF(I)) 

C2=(1-P)**NF(I) 
C3=C2/(1-P) 
C4=C3/(1-P) 
C5=C4/(1-P) 

C6=((1-P)^FTF(I))/P-FTNF(I) 
SUNKI1=SUMKI1+(C1*C1*C5)/P+NF(I)*NF(I)*C4 
SUMKI 2=SUMKI 2+C6*C6=*^P*C3+FTNF  ( I )  *FTNF  ( I  )  ^C2 
5      CONTINUE 

*  COMPUTE  STANDARD  DEVIATIONS 

PSD1=SQRT(1/SUMKI1) 
PSD2=SQRT(1/SUMKI2) 

*  COMPUTE  CONFIDENCE  LIMITS 

LL1=P-1.96*PSD1 
UL1=P+1.96*PSD1 
LL2=L-1.96*PSD2 
UL2=L+1.96'^PSD2 

*  OUTPUT 

WRITE  (6,320)  LL1,UL1 
320   FORMAT  ('0.95  ASSYMPTOTIC  CONFIDENCE  LIMITS  FOR  P:  ' ,F5 .4 , ' , ' ,F5 .4 
C) 
WRITE  (6,330)  LL2,UL2 
330   FORMAT  ('0.95  ASSYMPTOTIC  CONFIDENCE  LIMITS  FOR  LAMBDA:  ',F5.4, 
C,F5.4) 

RETURN 

END  ■    •  •    . 


I  I 
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APPENDIX  F 
QUANTILE-QUANTILE  PLOTS 
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Figure  F.l     [ixponential  Quantile-Quamile  Plot:  WUC  46X1600. 
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Figure  r.2     Lxponential  Quantile-Quantile  Plot:  WUC  4622100. 
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Figure  f'.3     Exponential  Quantile-Quantile  Plot:  WUC  632Z100. 
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Figure  F.4     Exponential  Qiiantile-Quantile  Plot:  WUC  74A1500. 
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Figure  F.6     IZxponential  Quantile-Quantilc  Plot:  WUC  5772200. 
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Figure  F.7     Lxponcntial  Quantilc-Quantile  Plot:  WUC  69 IS  100. 
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Figure  F.8     WeibuU  Quantile-Quantile  Plot:  WUC  46X1600. 
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Figure  F.9    WeibuU  Quantile-Quantile  Plot:  WLC  4622100. 
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Figure  F.IO    Weibull  Quantilc-Quantilc  Plot:  WL'C  632/100. 
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Figure  F.ll     WeibuU  Quantile-Quantile  Plot:  WUC  74A1500. 


86 


iXji  M     1-    i    i         \            \ 

C/) 

<•»/ 

\   •*V:        :           i"             i                i            :                           !                                     I 

VM* 

'      i       1  '^      '           :              1                '            !                          i                                 " ■ 

M  i\i    1     1     M         '           I 

o 
o 

!                     !                     '            '                           i\^'                                              i                                                    •                                      =                                                                                                                                                                                                        ^'^ 

1   M 1    !  \      1        1      i             i                \ 

M  M  i   \  1                1       ■   ^1 

•      i      ;    :       ;          •    ^L    '                                                                             *: 

M  Mj     j   \     i     i 

to 

^■+-i;i::t±:::::::^':::±=::;:::rz::7:r:-:f: 

—   _j 
< 

■t'1-Tt""T 1  "" '^ "\ i ''' *' 

Z 

■  1         III  y^.^-.^.— ....-.— .~. >......> T%t, ; _.« :..  If, 

Ul 

M    M   .j      i       1        i\j             i 

" 

"t                              K- 

^— 

L.I i....4...f..„_^___...i.„ «.„......._..> »\ „ __ ...„• 

UJ 

•       i             i          :              :                 ;                   ;              :*\                                                                ^ 

m 

;    i     i  1     ;        :         i           i        i  *\  '        '•■                      *\ 

u 

:       1       :    :         i            •               f                 :                      \ 

2 

.,      4. — p.,....-^           * ^ .f J ♦v.           ♦■ .• 

>— 

^*.   ...^,.   .f...  !..*....«... »*.._........^.M.  •...•.....!....    .....^.•..,...-  •- .  V^.... ....*.. ...f 

•       i        i    '         :             :                                  ^             '             *    \       ■ 

. 

i    i     ;  i      ;        1          i           •        ^               \*                   .    «: 

1 r"]"!""']-* ]■"■" i i i i\ " ; 

rH] \ i I \ j --V -; 

^Xp;::^^ 

d 

0>   0»    «nO     »f)        O           »D             O         ID                       •-                                  •- 

^  o»  c>cn    r^       iD        CM          —                                                     q 

o>                          S3iuNvno 

Figure  F.12    Weibull  Quantile-Quantilc  Plot:  WLC  74A5M0(.). 
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Figure  F.13     WeibuU  Quantile-Quantilc  IMot:  WLC  5772200. 
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Figure  F.14    WeibuU  Quantilc-Quantile  Plot:  WLC  69 IS  100. 


89 


APPENDIX  G 
DUAL  PARAMETER  OPTIMIZER 

1.        QUASI-NEWTON  METHOD 

The  Quasi-Newton  method  was  used  to  determine  the  maximum  likelihood 
estimates  for  the  WeibuU  and  Mixed  models.  This  method  was  derived  from  the 
Newton's  method  through  the  Taylor's  expansion  of  the  objective  function,  ({x),  where 
X  is  the  n  x  1  column  vector  containing  values  for  the  n  parameters  being  estimated. 
The  Newton  method  uses  the  quadratic  approximation  to  fl^x),  equation  G.l.  evaluated 
at  x' ^  to  iteratively  determine  an  optimal  solution. 

n;x)~fi:x(^Vvfl:x(^^)'^Ax  +  .5Ax'^V^fi;x(^W  (eqnG.l) 

The_  next  point  in  the  sequence  is  found  by  forcing  the  -gradient  of  the  quadratic 
approximation  (equation  G. 2)  to  zero  and  solving  for  x'         \    ^ 

VfIx)=VtIx(^b+V^fi;xW)Ax=0  '  (eqnG.2) 

Allowing  Ax  to  equal  x^       ^-x^  ^  equation  G.2  can  be  solved  for  x^         \ 

^{k+l)=^(k).y2f^^(k))-lyf(^(k))  (eqnG.3) 

The  two-parameter  optimizer  displayed  in  Appendix  G  did  not  use  the  Newton 
method  discussed  above  because  the  technique  can  be  unreliable  for  nonquadratic 
functions.    Instead,    a    more    robust    Quasi-Newton    method    was    used.    The    major 

difference  between  the  two  methods  is  that  the  Quasi-Newton  does  not  require  use  of 

7  1 

the  inverse  Hessian  matrix  V  f^x)    .    Instead,  it  uses  an  n  x  n  "A"  matrix  to  estimate 

the  inverse  Hessian  matrix.  The  A^^  matrix  was  iteratively  updated  using  the  Broyden- 

Fletcher-Shanno  method,  and  should  be  a  fair  approximation  to  V-f^x")    ,  when  x''  is 

the  optimal  solution.  Substituting  A^'^^  for  Vfi^x'   0      iri  equation  G.3  results  in  the 

Quasi-Newton  method  for  the  determination  of  the  next  point  in  the  sequence.    This 

iterative   process   continues  until   the   objective   function  is   maximized.     [Ref  12:  pp. 

112-114] 
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2.   FORTRAN  PROGRAM 

*  COMPUTES  MAXIMUM  LIKELIHOOD  ESTIMATES  FOR  TWO  PARAMETER 

*  MODELS  USING  QUASI-NEWTON  METHOD. 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/ RL/ ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 

COMMON/ INT/ CONFLG , COUNT , TFLAG , CENSOR , NF 

DIMENSION  CENSOR(200) ,NF(200) ,DIR(2) , FTNF(200) 

DIMENSION  FTF (200) ,GRAD(4) , PARAM(4) 

DIMENSION  A(4),A1(4),A2(4),DELP(2),DELG(2) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

*  INITIALIZE  'A'  -  APPROXIMATION  FOR  INVERSE  HESSIAN 

DATA  A(l)/1.0D0/A(2)/0.0D0/A(3)/0.0D0/A(4)/1.0D0/ITER/0/IRES/21/ 

CONFLG=0 

DIR(1)=-1.0D0 

DIR(2)=-1.0D0 

*  INPUT  DATA 

CALL  INPUT 

WRITE  (6,*)  ' COUNT ', COUNT 

*  OUTPUT  STARTING  VALUES  ; 

WRITE  (6,490) 

WRITE  (6,500)  (PARAM(I),I=1,2) 

*  COMPUTE  LOGLIKELIHOOD  VALUE  FOR  INITIAL  STARTING  ESTIMATES 

CALL  LLIKE(l) 

WRITE(6,510)  FV  .      ■ 

*  CHECK  FOR  CONVERGENCE 

CALL  CONCHK(l) 

IF  (CONFLG. EQ.l)  GO  TO  1000 

1      ITER=ITER+1 

WRITE  (6,520)  ITER 

*  IF  ITERATIONS  EQUAL  FIFTY,  PROGRAM  IS  TERMINATED  EARLY 

IF  (ITER. GE. 50)  THEN 

WRITE  (6,560) 

GO  TO  1000 
END  IF 

*  IF  IRES=1,  OBJECTIVE  FUNCTION  IS  NO  LONGER  IMPROVING  AND  GRADIENT 

*  IS  NOT  EQUAL  TO  ZERO.  THE  "A"  MATRIX  IS  RESET  TO  THE  IDENTITY  MATRIX. 

IF  ( ITER. EQ. IRES)  THEN 

A(1)=1.0D0 

A(2)=0.0D0 

A(3)=0.0D0 

A(4)=1.0D0 

IRES=IRES+20 
END  IF 

*  COMPUTE  DIRECTION  OF  IMPROVEMENT 

DIR(1)  =  -1'^(A(1)'^GRAD(1)+A(2)*GRAD(2)) 
DIR(2)=-1*(A(3)*GRAD(1)+A(4)*GRAD(2)) 

*  COMPUTE  LINE  STEP  SIZE  (ALPHA) 

CALL  LINE 

*  COMPUTE  DIFFERENCE  BETWEEN  PARAM(K+1)  AND  PARAM(K) ,  AND 

*  DIFFERENCE  BETWEEN  GRAD(K+1)  AND  GRAD(K) 

DO  10  1=1,2 

DELP ( I ) =PARAM ( 1+2 ) -PARAM ( I ) 
DELG ( I ) =GRAD ( 1+2 ) -GRAD ( I ) 
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10 


CONTINUE 


*  UPDATE  'A'  MATRIX 

C1=DELP ( 1 ) *DELG ( 1 ) +DELP ( 2 ) *DELG ( 2 ) 

Al ( 1 ) =DELP ( 1 ) *DELP ( 1 ) / CI 
Al ( 2 ) =DELP ( 1 ) *DELP ( 2 ) /CI 
Al ( 3 ) =DELP ( 1 ) *DELP ( 2 ) / CI 
A1(4)=DELP(2)*DELP(2)/C1 

C2=1-((DELP(1)'^DELG(1))/C1) 
C3=-1*((DELP(1)*DELG(2))/C1) 
C4=-1*((DELP(2)*DELG(1))/C1) 
C5=1-((DELP(2)*DELG(2))/C1) 

A2(1)=C2*(A(1)*C2+A(2)*C3)+C3'^(A(3)*C2+A(4)*C3' 
A2(2)=C2*(A(1)*C4+A(2)*C5)+C3'^(A(3)*C4+A(4)*C5' 
A2(3)=C4'^(A(1)*C2+A(2)'^C3)+C5*(A(3)*C2+A(4)*C3' 
A2(4)=C4'^(A(1)*C,4+A(2)*C5)+C5'*^(A(3)*C4+A(4)*C5; 

DO  20  1=1,4 

A(I)=A2(I)+A1(I) 
20     CONTINUE 


*  UPDATE  OLD  PARAMETER  ESTIMATES  AND  GRADIENTS  TO  NEW  PARAMETER 
■^   ESTIMATES  AND  GRADIENTS 


25 


DO  25  1=1.2 

PARAM (I ) =PARAM ( I +2 ) 
GRAD(I)=GRAD(I+2) 
CONTINUE 


*  OUTPUT  UPDATED  PARAMETERS 

WRITE  (6,500)  (PARAM(I),I=1,2) 

*  OUTPUT  FUNCTIONAL  VALUE 

FV=-1*FVAL 
WRITE(6,510)  FV 

*  CHECK  FOR  CONVERGENCE 

CALL  CONCHK(l) 

IF  (CONFLG.EQ.l)  GO  TO  1000 

*  START  NEW  ITERATION 

GO  TO  1 

*  PROGRAM  TERMINATION 
1000   WRITE  (6,530) 


WRITE  (6,540) 
WRITE  (6,550)  A 


*  OUTPUT  FORMAT 

490    FORMAT  ( ' OSTARTING  VALUES :',/,' 


) 


500 

FORMAT 

510 

FORMAT 

520 

FORMAT 

530 

FORMAT 

540 

FORMAT 

550 

FORMAT 

560 

FORMAT 

STOP 

END 

'OPHAT  IS:  ',F15.8,/',  '  LHAT  IS:  '  F15.8) 
'OLOGLIKELIHOOD  VALUE  IS:  ',F13.6) 
'OITERATION:  ',13,/,'  ') 

'OPROGRAM  COMPLETION') 

'OFINAL  "A"  MATRIX  IS : ' ) 

'0' ,2(F15.8),/,'  •,2(F15.8)) 

'OUNSUCCESSFUL  TERMINATION:  ITERATIONS  =50') 


SUBROUTINE  INPUT 
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*  READS  FLIGHT  AND  MAINTENANCE  RECORDS  USED  IN  THE  CALCULATION 

*  OF  MIXED  MODEL  MAXIMUM  LIKELIHOOD  ESTIMATES  FOR  PHAT  AND 

*  LAMBDA  HAT 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 

COMMON/ INT/ CONFLG , COUNT , TFLAG , CENSOR , NF 

DIMENSION  CENSOR(200) ,NF(200 ) ,DIR(2) , FTNF(200) 

DIMENSION  FTF(200) ,GRAD(4) ,PARAM(4) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

*  STARTING  POINTS  FOR  PARAM ( 1 ), PHAT ,  PARAM ( 2 ) , LHAT ,  AND  RECORD  COUNTER 

PARAM(1)=.67D0 
PARAM(2)=.01D0 
C0UNT=1 

*  INPUT  FLIGHT  DATA 

5     READ  (3,500,END=1000)  CENSOR ( COUNT ) ,NF( COUNT ), FTNF (COUNT ), FTF (COUN 
CT) 

C0UNT=C0UNT+1 
GO  TO  5 

*  INPUT  FORMAT 

500   FORMAT  (II , IX, 14, 5X,F6. 1 , 1X,F3 . 1) 

1000   C0UNT=C0UNT-1 
RETURN 
END 

SUBROUTINE  CONCHK(N) 

*  DETERMINES  IF  OPTIMAL  SOLUTION  HAS  BEEN  FOUND 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 

COMMON/ INT/CONFLG , COUNT , TFLAG , CENSOR , NF 

DIMENSION  CENSOR(200) ,NF(200) ,DIR( 2 ) , FTNF(200) 

DIMENSION  FTF(200) ,GRAD(4) ,PARAM(4) 

INTEGER  CENSOR, CONFLG, COUNT, TFLAG 

REALMS  NORMG 

*  TERMINATE  IF  THE  NORM  OF  THE  GRADIENT  (NORMG)  IS  LESS  THAN  THE 

*  THE  SPECIFIED  CRITICAL  VALUE  (CRIT) 

N0RMG=((GRAD(N)*GRAD(N))+(GRAD(N+1)*GRAD(N+1 )))**. 5 

IF  (NORMG. LE.l.)  THEN 

C0NFLG=1 

WRITE  (6,500)  ■ 
END  IF 

FTEST=ABS (UFVAL-FVAL) 

IF  (FTEST.LE. .000001. AND. NORMG. GT.l.)  THEN 

C0NFLG=1 

WRITE  (6,510) 
END  IF 


*  OUTPUT  FORMAT 

500    FORMAT  ( ' C 

510    FORMAT  ( ' OPROGRAM  TERMINATION  -  LACK  OF  FUNCTIONAL  IMPROVEMENT') 


599    E9B11M  r?!^?™L_. TERMINATION  -  OPTIMALFOUND^) 


RETURN 
END 
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SUBROUTINE  LINE 

*  COMPUTES  MAXIMUM  STEP  SIZE  (ALPHA) 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/ ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 

COMMON/ INT/ CONFLG , COUNT , TFLAG , CENSOR , NF 

DIMENSION  CENSOR(200) ,NF(200) , DIR(2 ) , FTNF(200) 

DIMENSION  FTF(200) ,GRAD(4) ,PARAM(4) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

REAL '^8  LB 

UFVAL=FVAL 
■  K=l 
TFLAG=0 
ALPHA=O.ODO 
ALPHAK=.0001D0 
LB=O.ODO 

*  COMPUTATION  OF  SEARCH  INTERVAL  TAU.  TAU  IS  APPROXIMATELY 

*  .618  AND  TAUl  IS  APPROXIMATELY  .382 

TAU=(5**.5-l)/2 

TAU1=1-TAU 

*  BOUNDING  PHASE  -  SWANN ' S  METHOD 

STEP=.00001D0 

PARAM(3)=PARAM(1)-STEP*DIR(1 
PARAM(4)=PARAM(2)-STEP^DIR(2 
CALL  LLIKE(3) 
F1=FVAL 

par;\m(3)=param(i)+step*dir(i) 

PARAM  ( 4 )  =PARAM  ( 2 )  +STEP'^DIR  ( 2 ) 

CALL  LLIKE(3) 

F2=FVAL 
IF  (F1.GE.F2)  then 

F1=F2 

GO  TO  5 
END  IF 
IF  (F2.GT.F1)  STEP=-1*STEP 

5     ALPHA=ALPHA+(2^^K)'^STEP 

PARAM ( 3 ) =PARAM ( 1 ) +ALPHA*DIR ( 1 
PARAM  ( 4 )  =P  ARAM  ( 2.)  +  ALPHA^D  I R  ( 2 

*  COMPUTE  LOGLIKELIHOOD  VALUE 

CALL  LLIKE(3) 
F2=FVAL 

*  IF  NEW  VALUE  IS  LESS  THAN  PREVIOUS  VALUE  TERMINATE  BOUNDING 

IF  (F2.GT.F1)  GO  TO  10 

F1=F2 

K=K+1 

GO  TO  5 

10     UB=ALPHA 

*  CONDUCT  GOLDEN  SECTION  SEARCH 

B1=TAU1*(UB-LB)+LB 

PARAM(3)=PARAM(1)+B1*DIR(1 

PARAM(4)=PARAM(2)+B1*DIR(2 

CALL  LLIKE(3) 

FB1=FVAL 

B2=TAU^(UB-LB)+LB 

PARAM(3)=PARAM(1)+B2*DIR(1 

PARAM(4)=PARAM(2)+B2*DIR(2 

CALL  LLIKE(3) 
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FB2=FVAL 
50    IF  (FB1.LT.FB2)  THEN 

ALPHA=B1 

UB=B2 

B2=B1 

FB2=FB1 

B1=TAU1*(UB-LB)+LB 

PARAM(3)=PARAM(1)+B1*DIR(1) 

PARAM(4)=PARAM(2)+B1*DIR(2) 

CALL  LLIKE(3) 

FB1=FVAL 

GO  TO  100 
END  IF 
IF  (FB2.LE.FB1)  THEN 

ALPHA=B2 

LB=B1 

B1=B2 

FB1=FB2 

B2=TAU*(UB-LB)+LB 

PARAM(3)=PARAM(1)+B2*DIR(1) 

PARAM(4)=PARAM(2)+B2*DIR(2) 

CALL  LLIKE(3) 

FB2=FVAL 

GO  TO  100 
END  IF 

*  CHECK  FOR  LINE  SEARCH  TERMINATION 
100    CALL  TERM(ALPHAK) 

ALPHAK=ALPHA 

IF  (TFLAG.EQ.l)  GO  TO  5000 

GO  TO  50 

5000   PARAM(3)=PARAM(1)+ALPHA*DIR(1) 
PARAM  ( 4 )  =PARAM  ( 2 )  +ALPHA'*^D  I R  ( 2  ) 
CALL  LLIKE(3) 

RETURN 
END 

SUBROUTINE  TERM(ALPHAK) 

*  DETERMINES  IF  CONVERGENCE  CRITERIA  IS  ESTABLISHED  FOR  LINE 

*  SEARCH  TERMINATION 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/ ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 

COMMON/ INT/ CONFLG , COUNT , TFLAG , CENSOR , NF 

DIMENSION  CENSOR(200) ,NF(200) ,DIR(2) ,FTNF(200) 

DIMENS ION  FTF ( 200 ) , GRAD ( 4 ) , PARAM ( 4 ) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

TFLAG=0 

CTEST=ABS ( (ALPHAK-ALPHA) /ALPHA) 

IF  (CTEST.LE. .003)  TFLAG=1 

RETURN 
END 

SUBROUTINE  LLIKE(N) 

*  COMPUTES  LOGLIKELIHOOD  FUNCTIONAL  VALUE  AND  GRADIENTS  FOR 

*  FOR  THE  WEIBULL  MODEL 

IMPLICIT  REALMS  (A-H,0-Z) 
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COMMON/RL/ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 
COMMON/ INT/ CONFLG , COUNT , TFLAG , CENSOR , NF 
DIMENSION  CENSOR(200) ,NF(200) ,DIR(2) , FTNF(200) 
DIMENSION  FTF(200) ,GRAD(4) ,PARAM(4) 
INTEGER  CENSOR, CONFLG, COUNT, TFLAG 

*  ESTABLISH  BOUNDARY  CONDITIONS 

*  PARAM (N)=SHAPE  PARAMETER; PARAM (N+1)=SCALE  PARAMTER 

IF  (PARAM(N) .LT. .0000000001)  PARAM(N)= . 0001 

IF  (PARAM(N+1) .LT. .0000000001)  PARAM(N+1 )= .0001 

FVAL=0 . 

GRAD(N)=0. 

GRAD(N+1)=0. 

*  INPUT  FLIGHT  DATA 

DO  5  1=1, COUNT 

C1=PARAM(N+1)*FTNF(I) 

C2=C1**PARAM(N) 

IF  (Cl.EQ.O.)  Cl=. 00001 

*  UNCENSORED  RECORDS 

IF  (CENSOR(I) .EQ.O)  THEN 

C3=PARAM ( N+l) * (FTNF ( I ) +FTF ( I ) ) 

C4=C3^*PARAM(N) 

C5=EXP(-C2' 

C6=EXP(-C4; 

*  FUNCTIONAL  VALUE  CALCULATIONS 

FVAL=FVAL-LOG(C5-C6)  ' 

*  GRADIENT  CALCULATIONS 

GRAD (N ) =GRAD (N ) + ( ( C5*C2*LOG( CI ) ) - ( C6*C4*L0G ( C3 ) ) ) / ( C5 -C6 ) 
GRAD  (N+l )  =GRAD  (N+l )  +  ( C5'^  ( FTNF  ( I )  '^^PARAM (N) )  -C6^  ( ( FTNF  ( I )  +FTF  ( I ) 
C    )'^^PARAM(N)))^((PARAM(N)*PARAM(N+1)*^(PARAM(N)-1)))/(C5-C6) 
GO  TO  5 
END  IF 

*  CENSORED  RECORDS 

*  FUNCTIONAL  VALUE  CALCULATIONS 

FVAL=FVAL+C2 

*  GRADIENT  CALCULATIONS 

GRAD (N ) =GRAD (N ) + ( C2*L0G ( CI ) ) 

GRAD(N+1)=PARAM(N)'^(FTNF(I)^*PARAM(N))*(PARAM(N+1)*^(PARAM(N)-1 
C  )) 

5      CONTINUE 

RETURN 
END 

SUBROUTINE  SD 

*  COMPUTES  THE  STANDARD  DEVIATION  FOR  THE  MAXIMUM  LIKELIHOOD 

*  ESTIMATORS   FOR  THE  WEIBULL  DISTRIBUTION,  WHICH,  IS  USED 

*  TO  CONSTRUCT  .95  CONFIDENCE  INTERVALS  FOR  THE  SHAPE  AND 

*  SCALE  PARAMETERS. 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON/RL/ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 

COMMON/ INT/ CONFLG , COUNT , TFLAG , CENSOR , NF 

DIMENSION  CENSOR(200) ,NF(200) ,DIR(2) ,FTNF(200) 

DIMENSION  FTF (200) , GRAD (4) , PARAM (4) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

REALMS  LL1,LL2 

SUMKI1=0. 
SUMKI2=0. 
C=(PARAM(1)*(PARAM(2)**(PARAM(1)-1)))'^'^2 
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INPUT  FLIGHT  DATA 
DO  5  1=1, COUNT 

IF  (FTF(I).EQ.O.)  FTF(I)=. 000001 

IF  (FTNF(I).EQ.O.)  FTNF(I)=-000001 

C1=(FTNF(I)+FTF(I)) 

C2=PARAM(2)'^C1 

C3=PARAM(2)*FTNF(I) 

C4=EXP(-1^(C2**PARAM(1)))/EXP(-1*(C3**PARAM(1))) 

SUMKI1=SUMKI  1+C*  ( ( ( ( FTNF  ( I )  **PARAM  { 1 )  -  CI  **PARAM  ( 1 )  )  '^C4 )  '*'  ( ( FTNF  ( I )  * 
C  *PARAM(1)-C1*'*^PARAM(1))*C4))/(1-C4)  +  (FTNF(I)**PARAM(1))*(FTNF(I 
C    )**PARAM(1))*EXP(-1^(C3**PARAM(1)))) 

SUMKI2=SUMKI2+((((C3*PARAM(1)'^L0G(C3))-(C2**PARAM(1)*L0G(C2)))^C4) 
C   **2)/(l-C4)+(C3**PARAM(l)*L0G(C3))*^2 
CONTINUE 

PSD1=SQRT(1/SUMKI1) 
PSD2=SQRT(1/SUMKI2) 

CONFIDENCE  INTERVALS 

•LL1=PARAM(1)-1.96*PSD1 
ULl=PARAM(l)+1.96*PSDl 
LL2=PARAM(2)-1.96*PSD1  • 
UL2=PARAM(2)+1.96^PSD1 

OUTPUT 

WRITE  (6,100)  PSDl,LLl,ULl 
WRITE  (6,110)  PSD2,LL2,UL2 


100   FORMAT 
110    FORMAT 

RETURN 
END 


('OPSDl:  ',F6.5,'  CONFIDENCE  INTERVAL  ( ' , F6 . 5 , ' , ' , F6 . 5 , ' ) ' ) 
('0PSD2:  ',F6.5,'  CONFIDENCE  INTERVAL  ( '  , F6 . 5 , '  , '  , F6 . 5  ,  ' ) '  ) 


SUBROUTINE  LLIKE(N) 

/^  COMPUTES  LOGLIKELIHOOD  FUNCTIONAL  VALUE  AND  GRADIENTS  FOR 

*  THE  MIXED  MODEL. 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/RL/ALPHA , FVAL , UFVAL , DIR , FTNF , FTF , GRAD , PARAM 

COMMON/ INT/ CONFLG , COUNT , TFLAG , CENSOR , NF 

DIMENSION  CENSOR(200) ,NF(200) , DIR(2) , FTNF(200) 

DIMENSION  FTF(200) ,GRAD(4) , PARAM(4) 

INTEGER  CENSOR , CONFLG , COUNT , TFLAG 

*  ESTABLISH  BOUNDARY  CONDITIONS 

IF  (PARAM(N).LT. .0000000001)  PARAM(N)= .0001 

IF  (PARAM(N) .GE.l.)  PARAM(N)= . 9999 

IF  (PARAM(N+1).LT. .0000000001)  PARAM(N+1 )= . 0001 

FVAL=0 . 

GRAD(N)=0. 

GRAD(N+1)=0. 

C1=1-PARAM(N) 

C2=L0G(C1) 

*  INPUT  FLIGHT  DATA 

DO  5  1=1, COUNT 

*  UNCENSORED  RECORDS 

IF  ( CENSOR ( I ).EQ.O)  THEN 

C3=EXP(-PARA>I(N+1)^FTF(I)) 

*  FUNCTIONAL  VALUE  CALCULATIONS 

FVAL=FVAL-LOG( 1-C3^C1 ) - (NF ( I ) -1 ) *C2+PARAM (N+1 ) *FTNF ( I ) 

*  GRADIENT  CALCULATIONS 

GRAD(N)=GRAD(N)-(C3/(1-C1^C3))+(NF(I)-1)/C1 
GRAD(N+1)=GRAD(N+1)-(FTF(I)'^C1*C3)/(1-C1*C3)  +  FTNF(I) 
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GO  TO  5 
END  IF 


*  CENSORED  RECORDS 

*  FUNCTIONAL  VALUE  CALCULATIONS 

FVAL=FVAL-NF ( I ) *C2+PARAM (N+1 ) *FTNF ( I ) 

*  GRADIENT  CALCULATIONS 

GRAD(N)=GRAD(N)+NF(I)/C1 
GRAD ( N+ 1 ) =GRAD (N+1 ) +FTNF ( I ) 
5      CONTINUE 


RETURN 
END 
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APPENDIX  H 
WEIBULL  SIMULATION  OF  COMPONENT  FAILURES 

*  VARIABLE  DEFINITIONS 

*  ALPHA :WEIBULL  SHAPE  PARAMETER 

*  AC:MATRIX  USED  TO  STORE  EACH  AIRCRAFTS  PREVIOUS  FLIGHT  HOURS 

*  WITHOUT  FAILURE  AND  UNCONDITIONAL  SURVIVAL  FUNCTION 

*  AS:AIRCRAFT  SELECTED  FOR  FLIGHT  EVOLUTION 

*  CPF: CONDITIONAL  PROBABILITY  OF  FAILURE 

*  DTIME: DEPLOYMENT  TIME 

*  FAILrMATRIX  USED  TO  STORE  SIMULATED  FAILURES  AT  SPECIFIED  TIME 

*  INTERVALS  FOR  EACH  DEPLOYMENT  SIMULATION 

*  FT: FLIGHT  TIME 

*  GFC: GENERATED  FAILURE  CRITERIA 

*  HM: HISTORICAL  MEAN 

*  HSD:HISTORICAL  STANDARD  DEVIATION 

*  LHATrMAXIMUM  LIKELIHOOD  ESTIMATE  FOR  WEIBULL  SCALE  PARAMETER 

*  NAC:NUMBER  OF  AIRCRAFT  USED  IN  SIMULATION 

*  NI:NUMBER  OF  TIME  INTERVALS  TO  BE  USED  IN  SIMULATION  (500  HR  INT.) 

*  NS:NUMBER  OF  SIMULATIONS  TO  BE  RUN 

*  NSEED1:SEED  USED  TO  DETERMINE  AIRCRAFT  CHOSEN  FOR  FLIGHT 

*  NSEED2:SEED  USED  TO  DETERMINE  FLIGHT  DURATION 

*  NSEED3:SEED  USED  TO  HDETERMINE  IF  A  COMPONENT  FAILURE  WAS  OBSERVED 

*  TGATE:TIME  GATE 

*  TIMEiACCUMULATED  FLIGHT  TIME 

*  UCPS: UNCONDITIONAL  PROFAILITY  OF  FAILURE 

DIMENSION  AC(70,2) ,FAIL(20,100) ,EST(20,4) 

INTEGER  AS 
REAL  *  4  LHAT 

*  INITIALIZE  •  . 

DATA  LHAT/. 0149/ ALPHA/. 8079/ 

DATA  NI/16/NS/100/NAC/48/HM/2.18/HSD/.9897/ 

DATA  NSEED1/749684/NSEED2/4683957/NSEED3/ 5872654/ 

*  COMMENCE  SIMULATION(I) 

DO  1  1=1, NS 
NF=0 
TIME=0. 
TGATE=500. 
DO  5  K=1,NAC 
AC(K,1)=0. 
AC(K,2)=1. 
5        CONTINUE 

*  SIMULATE  AND  RECORD  COMPONENT  FAILURES  USING  500  HR  INTERVALS 

DO  10  J=1,NI 

*  LRND  AND  LNORM  ARE  SUBROUTINES  IN  A  RANDOM  NUMBER  GENERATOR  PACKAGE 

*  CALLED  LLRANDOMII:  A  NON-IMSL  PROGRAM  AVAILABLE  AT  THE  NAVAL 

*  POSTGRADUATE  SCHOOL. 

*  GENERATE  UNIFORM  (0,1)  RANDOM  VARIABLE  TO  DETERMINE  WHICH  AIRCRAFT 

*  WILL  BE  FLOWN 

11  CALL  LRND (NSEED1,U, 1,2,0) 

AS=(U*(NAC-1))+1.5 

*  GENERATE  NORMAL  (0,1)  RANDOM  VARIABLE  TO  COMPUTE  FLIGHT  DURATION 
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CALL  LN0RM(NSEED2,RN, 1,2,0) 

FT=(RN*HSD)+HM 

IF  (FT.LT..5)  FT=.5 

IF  (FT.GT.6.)  FT=6. 

*  COMPUTE  THE  UNCONDITIONAL  PROBABILITY  OF  SURVIVAL  FOR  S+FT  HRS , 

*  WHERE  S  EQUALS  THE  TOTAL  FLIGHT  HOURS  FLOWN  SINCE  THE  LAST 

*  COMPONENT  FAILURE,  AND  FT  EQUALS  THE  FLIGHT  TIME  OF  THE  SORTIE 

*  BEING  FLOWN. 

UCPS=EXP(-1*((LHAT*(AC(AS,1)+FT))**ALPHA)) 
CPF=1-UCPS/AC(AS,2) 

*  GENERATE  UNIFORM  (0,1)  RANDOM  VARIABLE  TO  DETERMINE  IF  COMPONENT 

*  FAILURE  OCCURS 

CALL  LRND(NSEED3,GFC, 1,2,0) 

*  COMPONENT  FAILURE 

IF  (GFC.LT.CPF)  THEN 

NF=NF+1 

AC(AS,1)=0. 

AC(AS,2)=1. 

GO  TO  15 
END  IF 

*  COMPONENT  SURVIVES 

AC(AS,1)=AC(AS,1)+FT 
AC(AS,2)=UCPS 

*  UPDATE  TIME 

15  TIME=TIME+FT 

IF  (TIME.GE.TGATE)  THEN 

FAIL(J,I)=NF 

TGATE=TGATE+500. 

GO  TO  10 
END  IF 

*  NEXT  AIRCRAFT  SORTIE  TO  BE  SIMULATED 

GO  TO  11 

*  NEXT  TIME  INTERVAL 
10       CONTINUE 

*  NEXT  DEPLOYMENT 
1      CONTINUE 

*  CALCULATE  MEAN  AND  STANDARD  DEVIATION  FOR  THE  NUMBER  OF  FAILURES 
^    OBSERVED  AT  EACH  TIME  INTERVAL  (MEAN=EST( J , 1 ) ;STANDARD 

^    DEVIATI0N=EST(J,2)) 

DO  20  J=1,NI 

DO  25  1=1, NS 

EST(J,1)=EST(J,1)+FAIL(J,I) 
EST(J,2)=EST(J,2)+FAIL(J,I)*FAIL(J,I) 
25       CONTINUE 
20     CONTINUE 

DO  30  J=1,NI 
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EST(J,1)=EST(J,1)/NS 

EST(J,2)=((EST(J,2)-(NS*EST(J,1)*EST(J,1)))/(NS-1))**.5 
30    CONTINUE 

*  95  PERCENT  CONFIDENCE  INTERVAL  (LOWER  LIMIT=EST( J, 3) ; 

*  UPPER  LIMIT=EST(J,4) 

DO  35  J=1,NI 

HINT=1.96*EST(J,2) 
EST(J,3)=EST(J,1)-HINT 
EST(J,4)=EST(J,1)+HINT 
35    CONTINUE 

*  OUTPUT 

TIME=500. 

*  WRITE  (6,1000) 
DO  40  J=1,NI 

WRITE   (6,1100)  TIME,EST(J,1),EST(J,2),EST(J,3),EST(J,4) 
TIME=TIME+500.. 
40     CONTINUE 

*  OUTPUT  FORMAT 

*1000   FORMAT  ('1  TIME   MEAN    STD  DEV   LOWER    UPPER') 
1100   FORMAT  (F5.0,4(2X,F7.3)) 

STOP 
END 
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APPENDIX  I 
GRAPHICAL  COMPARISONS 


Figure  1. 1     Flight  Hour  Models:  Model  Set  WLC:734II1U0. 
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Figure  1.2     Geometric  Model:  Model  Set  \VUC:734II100. 
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Figure  1.3 


light  Hour  Models:  Validation  Set  \VUC:734H100. 
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Figure  1.4     Geometric  Model:  Validation  Set  \VUC:734H100. 


105 


• 

2 

o 

0 

\  ^ 

■ 

ii 

o 
.o 

'••  \  \ 

l.Ui 

^^ 

xw 

o 

*                                \ 

Q 

uo 

.o 

"■■■••  \  \ 

2 
U 

2?; 

o 
« 

\\\ 

Sggu 

>»^— 

o 

.  \  \V 

cii^:;^ 

§^- 

wcp 

i» 

q 

V      K 

c^ 

'J)  ■ 

\\ 

a  w 

"           \     * 

.-^ 

^    ■  '^-^ 

"•     \ 

o  c 
o  •« 

'■•••  \ 

"■■■■•••\\ 

a 
o 
a 

N 

O 

V\ 

.  f% 

1 — ■ — -^ — ' ' 

09                       0»                       OC                       02                       Of                        0 

snamivi 

Figure  1.5     Flight  Hour  Models:  Model  Set  WUC:40X1600. 
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Figure  1.6    Geometric  Model:  Model  Set  \VUC:46XI600. 
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Figure  1.7     Flight  Hour  Models:  Validation  Set  \VLC:-40X16OO. 
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Figure  1.8     Geometric  Model:  Validation  Set  WLC:46xl600. 
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Rgure  1.9     Flight  Hour  Models:  Model  Set  \VUC:4622100. 
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Figure  1. 10    Geometric  Model:  Model  Set  WLC:462210U. 
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Figure  1.11     Flight  Hour  Models:  Validation  Set  WLC:46221 
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Figure  1.12     Geometric  Model:  Validation  Set  WUC:4622100. 
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Figure  1.13     Flight  Hour  Models:  Model  Set  WUC:56X2500. 
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Figure  1.14    Geometric  Model:  Model  Set  WUC:56X2500. 
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Figure  1.15     Flieht  Hour  Models:  Validation  Set  \Vi;C:56X250(). 
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Figure  1.16    Geometric  Model:  Validation  Set  \VLC:56X2500. 
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Figure  1.17     Flight  Hour  Models:  Model  Set  \VUC:632Z100. 
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Figure  1.18     Geometric  Model:  Model  Set  WUC:632Zi<H). 
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Figure  1.19     Flight  Hour  Models:  Validation  Set  \VL"C:632Z100. 
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Figure  1.20    Geometric  Model:  Validation  Set  \VLC:632Z100. 
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Figure  1.21     Flight  Hour  Models:  Model  Set  \VUC:74A1500. 
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Fi'jurc  1.22     Geometric  Model:  Model  Set  \VUC:74A1500. 
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Figure  1.23     Flight  Hour  Models:  Validation  Set  \VLC:74A15O0. 
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Figure  1.24    Geometric  Model:  Validation  Set  \VLC:74AI500. 
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Figure  1.25     Geometric  Model:  Model  Set  \VL"C:74A5M00. 
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Figure  1.26     Flight  Hour  Models:  Validation  Sot  WLC:74A5M00. 
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Figure  1.27    Geometric  Model:  Validation  Set  WUC:74A5M00. 
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Figure  1.28     Flight  Hour  Models:  Model  Sot  \VLC:5772200. 
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Figure  1.29    Geometric  Model:  Model  Set  WUC:5772200. 
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Figure  1.30     Flight  Hour  Models:  Validation  Set  \VUC:57722UO. 
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Figure  1.31     Geometric  Model:  Validation  Set  WUC:577220(.). 
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Figure  1.32     Flight  Hour  Models:  Model  Set  WL'C:69 18100. 
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Figure  1.33     Geometric  Model:  Model  Set  WLC:6918100. 
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Figure  1.34     Flight  Hour  Models:  Validation  Set  WUC:6918100. 
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Figure  1.35     Geometric  Model:  Validation  Set  \VUC:691S100. 
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